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ABSTRACT 


A zonal Navier-Stokes model, developed by J.C. Wu, is 
installed and verified on the NASA Ames Cray X/MP-48 
computer and is used to calculate the flow field about a 
fek O02 airfoil oscillating in pitch. Surface pressure 
distributions and integrated lift, pitching moment, and drag 
Seetficients and integrated lift, pitching moment, and drag 
coefficients versus angle of attack are compared to existing 
experimental data for four cases and existing computational 
data for one case. These cases involve deep dynamic stall 
and fully detached flow at and below a freestream Mach 
number of .184. The flow field about the oscillating 
airfoil is investigated through the study of pressure, 
Pemetcity, local velocity and stream function. Finally, the 


effects of pitch rate on dynamic stall are investigated. 
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ie LNTRODUCTION 


Dynamic stall is a phenomenon that refers to an airfoil 
delaying stall beyond its static stall angle due to a rapid 
change in angle of attack. Associated with this is the 
generation of a strong vortex that appears at the leading 
edge of the airfoil which expands as it moves aft, causing 
large excursions in the pressure, pitching moment and lift 
the airfoil experiences (Figure 1, taken from [Ref. 1]). 

Because the dynamic stall angle generally occurs at much 
higher angles of attack than the static stall angle, the 
maximum lift the airfoil generates can also be much higher 
than for steady conditions. Unfortunately this is a 
transient condition, as the lift drops sharply when the 
vortex is shed from the trailing edge. However, if the 
mechanisms that govern the initiation and development of 
this vortex and the dynamic stall phenomenon can be 
understood and controlled, it could open the door to much 
more maneuverable and higher performing aircraft. Nature 
provides a prime example of what is possible in the 
dragonfly, which uses vortex energy recovery to help achieve 
its remarkable maneuverability. 

NASA currently has a flight test program utilizing 
vortex generating leading edge slats to gain further insight 


into this fascinating area. Parallel studies are underway 
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by numerous groups to investigate dynamic stall via 
computational methods. The approaches that have received 
the most attention have used Navier-Stokes modeling, which 
has the flexibility to describe many flows. 

Navier-Stokes modeling has usually been formulated along 
velocity/pressure lines. This formulation has two major 
problems [Ref. 2]. The first problem is caused by the size 
of the flow field, which in mathematical terms, is infinite 
in extent. Boundary conditions at infinity must be 
Satisfied and an infinite flow field has to be modeled by a 
finite number of grid points. To do this, the boundary 
conditions must be previously known and multiple boundary 
condition solutions computed and compared to the _ known 
solution. Alternatively, a coordinate transformation can be 
made from the infinite flow field to a finite one. 
Unfortunately, this introduces a requirement for additional 
computations of an increasingly complex nature. 

The second problem is the large number of data points 
required to adequately model the flow field. Even for two 
dimensional cases this can be a significant impediment. 
Various grid spacing methods have been used to concentrate 
data points in the regions around the airfoil where high 
flow variable gradients occur, but the number of data points 
required is still quite high. 

An alternative method is to use velocity/vorticity 


modeling. By taking advantage of certain features of 


velocity and vorticity fields for incompressible viscous 


flow, 


Jee 


These 


certain conclusions may be reached [Ref. 3]. 


Vorticity can be neither created nor destroyed in the 
IntemLoirs OE thom me ruaic 


The total vorticity in the infinite unlimited space 
jointly occupied by the fluid and the solid is always 
ZeEEO. 


The rate of convective transport of Vorterera ae 
Eintce. The rate of diffusive transport is 
effectively finite. 


At large distances from the body, the velocity field 
approaches zero with increasing distance from the 
solid. 


At large distances from the solid, the vorticity field 
decays exponentially with increasing distance from the 
Sie) laisle 


The velocity field is uniquely determined by the 
vorticity distributed in the infinite unlimited space 
jointly occupied by the fluid and the solid. 
Alternatively, the velocity field is uniquely 
determined by the vorticity distribution in the fluid 
and the velocity condition on the solid boundary. 


features have three major salutary effects. [Ref. 4] 


The actual number of grid points that need to be 
solved for will generally be much less than the total 
number of grid points needed to define the flow field. 
This is because vorticity is generated solely at the 
solid boundary, and this vorticity is diffused only a 
short distance from the solid before being carried 
away by convection and diffusion. Therefore, the flow 
will be inviscid a short distance ahead of the solid 
and often inviscid at small to moderate distances 
above and below the solid. 


By taking advantage of the integral representation 
that the vorticity/velocity equations can be expressed 
in, the various regions of the flow field can be 
treated separately, with separate grids and distinct 
computational procedures that consider the length 


scales and defining equations, with no loss of 
accuracy and without a need to match the solutions of 
the various zones. Ultimately, an elegant method to 


compute flow field solutions becomes available that 


can accurately compute dynamic stall conditions 
without resorting to brute force velocity/pressure 
finite difference computations. 

3. The solution may be expressed as an integral which can 
be converted to a Fourier series expansion. [Ref. 5] 
This allows for very accurate computational solutions 
at each grid point. 

Currently the model is limited to two-dimensional 
incompressible flows, but the concepts are applicable to 
three-dimensional and compressible flows as well. 

A zonal Navier-Stokes model has been developed by J.C. 
Wu and his associates at the Georgia Institute of Technology 
and was made available for this’ study. An explicit, 
integro-differential methodology procedure is used to solve 
two-dimensional, Reynolds averaged, incompressible Navier- 
Stokes equations. 

The goals of the present study were: 


1. Install and verify the code on the Cray X/MP-48. 


2. Compare the code's solutions with previous 
experimental and theoretical results. 


3. Modify this code to enhance its utility. 


Il. MATHEMATICAL FORMULATION 


A. GOVERNING EQUATIONS DEVELOPMENT 
In the study of fluia@ yeiows, normally certain 
assumptions will be made. These assumptions allow emphasis 


on the specific features of interest, while providing useful 


simplifications in the theoretical development. Properly 
chosen, they can also promote the ease of numerical 
formulation and consequent solution. Accordingly, the 


assumptions made for the flow type in this study and their 
respective purposes are: 


1. Two-dimensional. The primary benefit of this 
assumption 1S a very significant reduction in 
computational requirements. A relatively minor amount 
of information is lost for spanwise flow and tip 
effects, with the salient aspects to be investigated 
remaining intact. 


2. Incompressible. A useful simplification to highlight 
the dynamic stall effects and allow for a more readily 
formulated computational model. 


3.  Viseouse A necessary element to describe the solid 
body and flow field interaction. This also helps to 
control the extent of flow field computations 
required. 


4. Unsteady. Required to allow for airfoil pitching. 


5. Turbulent. Required to describe the flow field under 
the conditions of interest. 


6. Reynolds averaged. A formulation of the Navier-Stokes 
equations for turbulent conditions. 


1. The Navier-Stokes Equations 


For three-dimensional unsteady flow in the x- 


direction: 
Du _ xy - iL oP, yu 
Pape 7 Pe 0 Ox Bi 
du du du du 
a = pein Rsiioad eit 
Oo p (ae + x + VY 3y* re 
2 2 2 
=px- 2 By BPR AE oS 
© ery oy OZ 
where: 
Du _ 9u ou ou ou _ £ 
at = aa De +v xy + W a5 substantial derivative of u 
ox = body forces 
op 
‘x = pressure forces 


5° su 3-u 
wage + a) + —l = viscous term. 
P ox dy OZ 


2. Reynolds Averaging for Turbulent Flows 
The boundary layer equation of motion in two 


dimensions is [Ref. 6]: 
» e) du 
9 = - SE + wale Sy) 


letting the following variables be defined as: 


where: 


Then: 


where: 


U(x,y) + u'(x,y,?t) 


Cc 
I 


V(x,y) + v'(%,y,t) 


< 
li 


B= Fy) ee a 


U(x,y) = mean x-direction velocity value 
over time 


u' = fluctuation Value ore velocity im 
x-direction 


i ace 


p = pressure 


Substituting u, v and p into the boundary layer 
equation, 


Taking the mean value of each tern, 


Observing that the mean value linear terms in a 
fluctuation component vanishes, 


Assuming that the derivatives of the mean value linear 
term vanishes, and 


Neglecting turbulent normal stress compared to the 
shearing stress, 


JU dU CE Omiya 
0 (U aaa V ay) a y | ou) 
-ou'v' = the Reynolds stress. 


Generally, the neglected terms are much smaller than -pu'v', 
and it should be noted that u' and v' can reach values as 
hae aS «6.1 «UZ: However, for the flow under consideration, 
these are still valid approximations. 
3. Prandtl's Mixing Length 

From the Reynolds stress development, another 
concept can be presented [Refs. 7,8]. Boussinesq introduced 
a mixing coefficient that related the Reynolds stress to the 
derivative for U. This is possible because for there to be 
a net momentum transfer from the higher to lower momentum 
layers -pu'v! PrOwonL ah octher words, uu’ and v' must be 
positively correlated. This is done by setting the Reynolds 


stress equal to the derivative of U or: 


where: 


Vv, = eddy viscosity which is a turbulent mixing 
Goefficrent. 


Peace moOoSetvedmrnaueydepends on U and is not a 
mioperty of the fluid. A relation between Vv, and the mean 
velocity was needed. Prandtl's idea was that if small 
"lumps" of fluid move from a lower to a higher average 


velocity location, the difference in its velocity compared 


to the surrounding mean velocity could be given by: 


where: 


G 
li 


f(y) for the simplified case, 


ro 
il 


Prandtl's mixing length. 


The physical significance of % is the distance in the y 
Girection the small lump of fluid must travel so that the 
difference between its velocity and U is equal to v; which 
is the y-direction fluctuation velocity at that location. 


After further development, Prandtl observed that: 
-opu'vi = pk (3y) = turbulent shearing stress 
This relation is very useful in the calculation of turbulent 


flows. 


B. GOVERNING EQUATIONS 


The equations of motion for the flow field under 


consideration are: 


Velocity > CHOU Cl 

field: Vev = 0 or yx + ay @) 

VOriEtCatey, ; > > ™ (ov 2 uy = ‘ 

frela: A eae oO Ox y 

Navier-Stokes 2) 2 
in the ou 4, CLierer so = = = <P 4 ee gu 
x-direction: gt OX Y p aX Jy 


10 


which, after taking the curl of both sides of the equation 
and applying the definitions of continuity and vorticity 


becomes [{Ref. 9]: 


OW (weV)v -— (veV)w + ae 
dt 
where: 
(w°V)v = stretching and rotation 
(v-V)w = convection 
v4 = diffusion 


It should be noted that for turbulent flows, Vv should be 


Beplaced by vo, OF: 


where: 


<— 
I 


ls effective viscosity 


= 
ll 


_ eddy viscosity 


H= kinematic viscosity. 
e 


< 
ll 


C. ANALYTICAL FORMULATION 
1. Kinematics 
Kinematics is the branch of dynamics which deals 
with the motion of bodies without reference to the forces 


acting on the bodies. Of the three governing equations, the 


ia 


two that comprise the kinematics of the flow are the 
equations of continuity andivewercie Together, they 
express the relationship between velocity and vorticity 
throughout the fluid at any point in time. 

There are two noteworthy aspects of these relations 
f[Refs. 10,11). 


1. The differential equations are linear, hence they can 
readily be formulated for solution by computational 


methods. 
2. The stress-strain relation is not a factor in@ene ce 
equations. This allows the fluid and solid to be 


treated together as one kinematic system which greatly 
Simplifies the computations. 


Through the use of fundamental solutions [Ref. 10] 
the continuity and vorticity field equations can be 


expressed as [Ref. 5]: 


>> > > > > > 
Wie =) ea J % xVoPdR, + $ ((vq"M9) = (Vo *Ny) J xV PCB, 
where: 


—_> 
the subscript "0" refers to the rg space, 
== 


Yq Space is defined by the vorticity field, OF, Wy= W (Ty, 


B is the boundary of the region R, 


> 
n is the outward unit normal vector on B, and 


di 1 
P= - 5, ln Eo which is the fundamental 
ro 
solution for a two-dimensional Poisson 
equation. 


The consequences of this formulation include: 
1. The first integral will be zero for the inviscid 


region, therefore it needs to be computed only in the 
viscous region. This greatly reduces the 


eZ 


computational time required, as the viscous region is 
Cwe woemmaaction Of the total flow field. 


The integrals can be expanded by Fourier series. This 
will provide more accurate solutions at each grid 
point than is possible via a straight finite 
difference method. 


The attached viscous and detached viscous zone 
solutions may be computed separately. This allows an 
optimized grid spacing to be employed in each Zone, 
where the respective length scales vary greatly. 
Figure 2 presents a flow zone portrayal. 
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Figure 2. Flow Zones 


Repeated computations and comparisons with either the 
boundary values or between the different flow zones 
are not required. 


The question of modeling a strong viscous-inviscid 
interaction is obviated completely. 


It should be noted that all of these aspects either 


reduce the number of calculations required, or enhance the 


ecuracyeor the solution. The net result is an elegant 


i, 


solution to the conundrum encountered by those who wish to 
model this type of flow. 
2 7 obec aes) 

Kinetics is the branch of dynamics which deals with 
the effects of forces on the motion of bodies. The 
vorticity transport equation falls into the realm of 
kinetics. Tt is nonlinear and elliptic in space. This 
requires knowledge of the solution for the entire boundary 
Conditions. For the outer boundary, this is met at the 
surface just inside the inviscid region. This boundary can 
change with each time step, but from the initial conditions, 
the vorticity will be zero along this surface. The interior 
boundary is the solid, and the vorticity values on this 
surface will need to be computed each time step. Because 
the vorticity values on the surface of the solid are not 
independent of the vorticity values in the interior of the 
fluid, this will have to be solved iteratively. 

3. Grid Generation 

In computational fluid dynamics, grid Cchoilcomiere 
something of an art. There are a number of conflicting 
goals to be considered. These include: 


1. An adequate number of grid points to accurately 
represent the flow conditions. 


2. A small enough number of grid points to help moderate 
the computational requirements. 


3. Proper grid spacing demands including a fine mesh in 
the regions of high flow variable gradients and a 
coarser mesh in the other regions. The latter will 


help moderate the computational resources required. 
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4. A grid that is easy to generate. 
5. A grid that simplifies and speeds computations. 

Fortunately, due to the present mathematical 
formulation, a very efficient grid scheme may be employed. 

By choosing the computational grid to be circular 
with radial lines, grid generation and computation goals can 
be met. Because the flow zones can be computed separately, 
grid spacing can be optimized for each Zone. Due to the 
relatively small distances from the solid that the vortical 
flow is transported to, the size of the grid required can be 
moderated. The only drawback to this method is that the 
computational plane must be conformally mapped onto the 
physical plane. This limits the airfoil selection to those 
that can be accurately mapped between planes. Fortunately, 
a number of airfoils are available via Joukowski transforms, 
including the NACA 0012, for which there is a wealth of 
previous data available for comparison. 

A salutary effect of using a Joukowski transform 
between the computational and physical planes is that the 
grid density in the physical plane is concentrated radially 
about the leading and trailing edges while becoming more 
Sparse radially over the upper and lower surfaces of the 
Spec EoOll . 

Lastly, the computational grid is body fixed, which 
eliminates having to recalculate the grid for each time 


step. Representative grids are included in Figures 3-6. 


1s) 





Figure 3. Boundary Layer Grid 
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4. Inifial and Boundary (Condieions 
Tnitially, the solid and the fluid are at rest iia 


solid is then impulsively started, and the discontinuity 
between the solid and the fluid generates a vorticity sheet 
at the boundary. As time progresses, vorticity is diffused 
away from the solid and transported into the fluid by both 
diffusion and convection. 

The four boundary conditions that must be met are 
for velocity and vorticity at the outer surface onmiene 
viscous region and at the solid/fluid boundary. As was 
previously mentioned, the outer vorticity boundary condition 
is met by setting the boundary just inside the inviscid 
region. The outer velocity boundary condition will simply 
be U. 

The inner VGrGtenty boundary Condi trons as 
previously stated, will need to be computed iteratively each 
time step, with the inner velocity boundary condition found 
by first determining the tangential velocity component, 
which is due to the solid's speed, angle of attack and 
oscillatory motion [Ref. 5} and then using that value and 
the relation v = Veoliq(rrt) (Ref. 13], to compute the inner 


velocity boundary "condiezone 
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Tl ioe rnoCREPTLONSORSTHE CODE 


The entire runstream has been named ZETA, which is an 
acronym for zonal procedure for evaluating turbulent and 
laminar flows. ZETA is composed of three primary sections: 

1. GEOM--grid generation and transformation, 
2. ZONST--main program in ZETA runstream, and 


3. Plotting routines--consisting of pressure computations 
ance various plotting options. 


miescomputcational loop in ZONST consists of [{Ref. 5]: 

1. Computing interior vorticity values by using the 
vorticity and velocity values from the previous time 
step in the vorticity transport equation. 

2. Computing new boundary vorticity values. 

3. Computing new velocity values. 

Appendix B contains the version of the Wu code used for 
this study and Appendix C contains notes on its employment 
at NASA Ames. 

A complete cycle through 40 degrees of pitch change at a 
reduced frequency of .15 requires 950 time steps and 
approximately 800 seconds of CPU time on the Cray X-MP/48 
Seiad average of .842 seconds per time step. ~ Initializa- 


tion can be completed with 25 time steps. These are much 


more modest requirements than codes of Similar ability. 


A. GEOM 
GEOM begins by reading the transformation and grid 


parameters. The transformation parameters define the 


ig 


airfoil that will be used while the grid parameters define 
the size and stretching coefficients which controler 
radial grid spacing. 

GEOM next generates the computational plane and then 
conformally maps it onto the physical plane. The outputs 
from GEOM are used by all the other programs. Grid pau 


spacing may be checked here and adjusted as necessary. 


B. ZONST 

ZONST is the main) program) of we4e te It uses the 
governing equations to compute the vorticity and velocity 
fields and generates the output used by the plotting 
rowednes. 


It starts by reading the grid information from GEOM and 


its own input parameters. The input options available 
inelude: 
Li eA rior temo erton This can be defined as a constant 
angle of attack, a rapidly pitching or an osci/lacwa 
motion. All pitching airfoil flow solutions require 


data from the appropriate constant angle of attack 
steady state case to be stored as part of the initial 
conditions. 


2. Time specifications allow the user to control the 
number of time steps to be run; the increment of the 
time steps and the time step values for which 
numerical and graphical output will be generated. 


3. Flow zone specifications can be used to control the 
boundary layer region where Navier-Stokes computations 
will be used. Generally, there will be no increase in 
accuracy but a significant increase in computation 
time when using Navier-Stokes vice boundary layer 
calculations in the boundary layer region. 
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4. Under-relaxation parameters allow tailoring of the 
computational methods to help optimize computing 
efficiency for various flow conditions. 

ZONST can be run for a specified period, checked and 
restarted. Should divergence occur, the last computed time 
step will be printed for trouble shooting. 


The output from ZONST is numerical, and via plotting 


routines, graphical. 


fe PLOTTING ROUTINES 

Due to the primary variables used in ZONST, the 
aerodynamic loads caused by the pressure distribution are 
not directly available. LOADS uses the vorticity 
information and Fourier series expansions of the vorticity 
integral representation to derive values for coefficients of 
pressure, lift, drag and moment for each angle of attack. 

The first three plotting routines listed use DISSPLA for 
their operating software, while the last one uses PLOT3D. 

Plotl portrays the physical plane vorticity grid, the 
boundary layer grid, flow zone demarcations and wake and 
turbulence grids. 

Plot2 generates streamline and vorticity contours as 
well as numerical output of the contours and grid points 
crossed. 

Pl@es. 15 the Hoads plotting progran. CoefLEici1ent of 
pressure versus non-dimensionalized chord length for each 


selected angle of attack may be displayed, or plots for the 


Zi 


coefficients of lift, drag and moment versus angle of atta 


for unsteady cases are available. 


Ze 


iso GLO LANDsDISCUSSION 


A. COMPARISON WITH EXPERIMENTAL DATA 

The first comparison was with experimental data from 
exert. 14) for a series of three different reduced 
frequencies, at the same Mach and Reynolds numbers. This 
was selected to highlight the time history dependent nature 
of dynamic stall. The plots pertaining to this comparison 
are enclosed in Appendix A. Considered were Cp versus non- 
@amensionalized chordlength or x/c, as well as Cy, Cy and Cp 
versus incidence angle, a. 

Reduced frequency is a pitching rate parameter and is 


defined as: 


where: 
W = circular frequency, 
C = airfoil chord length, and 
U = freestream velocity. 
The experimental data were taken at M = .072 to provide 


a valid approximation of incompressible flow. 
aed Seo? sche plots Of C,, Cp and €, versus &% pro- 


vide a good overview of the conditions the airfoil 


Zo 


experiences. In all three plots, the theoretical data 
follows the trends of the experimental data well, while the 
relative magnitudes and short term stall recovery features 
are less well represented. 

The theoretical values of Cy also track relatively well 
with the experimental data, though slightly under-represent= 
ing the pressures recorded experimentally. The theoretical 
results indicate the airfoil as beginning to stall at a 
slightly lower angle of attack than the experimental results 
do. At © = 20.8, the theoretical data show the vortex 
bubble as having been shed from the leading edge of the 
airfoil. Its progression across the upper surface can be 
followed in subsequent plots. 

Streamline and vorticity contour plots are presented to 
illustrate the physical phenomena but are not directly 
compared with experimental data. Vorticity buble 
generation and propagation are observed in each of the three 
cases. 

For the plots of rf = .15, stall onset is noticabie 
delayed from the rf = .099 case. Stall now occurs atvon 
24° for both the experimental and theoretical cases. As in 
the rf = .099 case, the trends show good correlation with 
post stall effects being less well modeled. 

For rf = .248, the plots of Cy, Cp and C,, versus @ jhe 


somewhat poorer correlation than the previous two cases. 
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There is greater offset between data, but stall occurs at 
a = 25° for both data sets on the Cy;, versus curves. 

For these three cases, there is a slight trend towards 
higher maximum values on the Cy; versus % curves with higher 
reduced frequency. Significantly, the theoretical results 
accurately reflect the trend towards a higher’ stall 


incidence angle with higher reduced frequency. 


B. COMPARISON WITH OTHER THEORETICAL DATA 

For the last case, a comparison was made between 
experimental data and two types of theoretical data. 
Experimental data are from [Ref. 14] and one version of the 
theoretical data is from (Ref. 15]. The conditions for this 
case were: reduced frequency = .199, mach number = .184 and 
oes 0 < 18°. 

Figure 7 is the experimental and other theoretical data. 
Figure 8 is the theoretical data from the Wu code. 
Referencing the Cp versus x/c plots in Figs. 7a,b and 8c, 
the theoretical data from the Wu code demonstrate a very 
high degree of accuracy with the experimental data, while 
the theoretical data from the other source show a decidedly 
less accurate picture. The other source indicates dynamic 
stall and post stall conditions when experimentally, the 
airfoil never stalled. 

Similar and pronounced differences also occur in the Cry, 
Sicha G. versus o& plots. As in the previous series, the Wu 


code data slightly under-represents the amplitude of the 
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pressure but still follows the trends well, while the other 
theoretical data which incorrectly indicates a dynamic stall 
eonaition. 

For this case, the data from the Wu code data are a more 
accurate reflection of the experimental data than the data 


from the other theoretical source. 


fo ENHANCEMENTS TO THE CODE 

There were two primary enhancements that were added to 
the Wu code. These are: 

1. The ability to output plotting data at non-regular 
intervals, to add flexibility to the code and to 
Simplify and facilitate comparison with experimental 
data. 

2. The availability of the code to generate physical 
plane velocity data to allow the use of Plot3D and 
its associated graphics functions. 

The modifications to provide physical plane velocity 
information centered around the transformation matrix that 
is used in GEOM to obtain the physical plane grid from the 
computational plane grid after the computational grid is 


generated. 


Tne scale factor of the transformation is defined as 


where Z is the physical plane and 


24 = relé 


£(7c) 


px) 


Z = £(oet?) 


while c is the computational plane, where the airfoil is 
represented by the unit circle and cylindrical coordinates 


are employed. 


Then, denoting = as DZDW, HSTAR is defined as: 
H* = complex(real(DZDW) ,-absolute imaginary (DZDW) ) 


HSTAR is written onto tape 20 and passed to ZONST, where 
it is used in the subroutine KMTCS. In KMTCS, ewe 
additional arrays are defined, VTOTCO and VTOTPH, which are 
the total computational plane and physical plane velocities, 


respectively. They are defined as: 


VEOTES 


complex (W1,W2) and 


VIOTEre 


VTOTCO/HSTAR 


where Wl and W2 are the computational plane velocities in 
the rotating frame of reference. Additional minor modifica- 
tions to the code were made to implement these changes. 
Representative plots are presented for four cases where 
the reduced frequency is .150 and Reynolds number is 1 10°. 


1. Laminar flow with the airfoil at a steady state 5.0° 
a. This is portrayed in Fig. 9. 
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Figure 9 
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2. Flow conditions immediately prior to stall onset are 
shown in Fig. 10. 


3. Initial indication of reverse flow at 21°7506° © ines 
ifvaleee This appears on the upper surface near the 
leading edge as the vorticity bubble starts to forn, 
which quickly spreads over the entire upper surface of 
the airfoil. 

4. Turbulent flow with the aingfoil at 23.890, Shore 
after the onset of dynamic stall. This is portrayed 
in Fig. 12. In all cases, alternate radial grid lines 
were deleted from the plots to enhance clarity, while 
a Similar display progression is followed for the 
first and last cases. 

Two additional modifications would help interpretation 
of the information. First, computing VTOTPH for the entire 
physical plane, not just in the viscous region. Second, 
rotating the grid so that the airfoil is graphically shown 
at its true angle of incidence. Both of these deficiencies 
are presented in Fig. 9a, and neither affects the accuracy 
of the graphical output in the viscous region. 

In Fig. 9a, the primary portion of the viscous Teq@en 
and part of the inviscid region are shown. 

Subsequent portions of Fig. 9 show increasing resolution 
of various portions of the plot im Fig.) oa- To aidwan 
comparison, the horizontal scale is consistent throughout 
the velocity plots. 

Fig. 9c shows the leading edge stagnation point and 
tangency of the boundary layer at the surface of the 


ALE ro] ie Fig. 9d is the midchord upper surface region. 


Fig. 9e is the aft portion of the airfoil. Fig. 9f is a 
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Figure 12 


detail plot of the boundary layer near the mid-span, upper 
surface. 

Fig. 11 1s the case for dynamic stall. Flow reversal is 
apparent over the entire upper surface, with the vortex 
bubble center indicated by the zero velocity vector located 
above approximately .3 chord. In Fig. lic, the @igmgag 
velocity gradients are readily apparent, from the leading 
edge around to the upper surface. 


More complex flow patterns can be readily portrayed. 
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V. CONCLUDING REMARKS 


The Wu code holds much promise to help unlock some of 


the mysteries of dynamic stall. Its strengths include the 
following: 

- Enough speed and efficiency that it can be operated ona 
VAX or similarly-sized computer. 

- Good success at indicating the trends of the cases 
studied. 

- Relatively accurate results. 

- Powerful diagnostic tool which can become a predictive 
cool . 

- Not Reynolds number limited. 

Other aspects that should be noted are: 

- The mathematical formulation iS more involved than a 
straight finite differencing of pressure/velocity 
Navier-Stokes equations. 

- As currently formulated, it is not readily applicable to 
arbitrary geometries, but a useful selection of 
Joukowski transforms is available. 

With the addition of compressibility and transition 
modeling and _ later, extension to three dimensional 


representation, its utility will continue to be enhanced and 


ts 


applications expand. 
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APPENDIX A 


REDUCED FREQUENCY VARIATION PLOTS 


Appendix A contains plots for variation of reduced 
frequency. For each of the three reduced frequencies, the 
order is: Cy, Cp and Cy, versus o 

Cp versus x/c with varying «o 
Streamline and vorticity contours with varying a. 
Theoretical: 


Experimental: -~------------ 


42 


pUstolJJeo0) YW 





yUsTOYya0) deiq 


Oo 
o) © >) 


=10).3 
—0.4 


yUaTOI{Jao0y) JUSUIOW 





-0.5 


10.0 


9.0 


0.0 


43 


0.4 


0.5 


X/C 


44 


a 
t 
rf 
Re 


=_ 
= 
— 
—_ 
— 
— 
— 
—_ 








45 





46 





00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 08 0.9 1.0 


47 


a 
t 





48 


aX 
ie 





0.0 0.1 0.2 O33 0.4 0.5 0.6 0.7 0.8 0.9 10 
X/C 


t 
rf 
Re =1.0*10° 





49 


OX 22.9 

t 69.520 
tet 0.099 

Re =1.0*10° 





50 


20.0 


15.0 





0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
X/C 


14.5 
t 87.120 
rf = 0.099 
Re =1.0*10° 





(omeol "O25 0S O44 05 OG O07 08 O09 10 
X/C 


aL 


ro | 

t 

rf 

Re =1.0*10— 





Diz 





—Cp 


20.0 


15.0 


10.0 


0.0 


25 


OX 

t 

Teli 

Re =1.0*10° 





Streamlines 


O 5.239 
t 4.000 
i 0.099 
Re =1.0*10° | 


Vorticity Contours 





54 


Streamlines 


Vorticity Contours 


OX 0.946 
t 8.000 
ayy 0.099 
Re =1.0*10° 





Sy) 


Streamlines 


Cx . 

8 
| rf 0.099 
ea Ogee 





56 


Streamlines 


Vorticity Contours 


X 

t | 

rf 0.099 
Re =1.0*10° 





51) 


Streamlines 


‘Vorticity Contours” 


10.432 





Re =10+10" | 


58 


Streamlines 


Vorticity Contours 





Se, 


Streamlines 


OX : 
t 
Ie 





60 


Streamlines 


Vorticity Contours 





61 


Streamlines 


O 

t 

et 

Re =1.0*10° 





62 


Streamlines 


Vorticity Contours 


et 20.827 
t 40.000 
rf 0.099 
Re =1.0*10° 





63 


Streamlines 


a a 


~ 
t 


roti , 
Re =10710 4] 





64 


Streamlines 


‘at 
1 
Ak 





65 


Streamlines 


a 
t 


O 
IU 





66 


Streamlines 


Vorticity Contours 





67 


Streamlines 


OX 


6 60.000 
fol 0.099 
ke =1.0°10_ 





68 


Streamlines 





69 


Streamlines 





10, 


Streamlines 


Vorticity Contours 





7A 


Streamlines 





72 


streamlines 


Vorticity Contours 





oe) 


Streamlines 


168072 

84.000 
ie 0.099 
= LOMO” 





74 


Streamlines 


Vorticity Contours 





Iae: 


Streamlines 


—- 
— 
——_ 
—_—= 
=-—— 
—«z 
== 
— 


-Vorticity Contours ~ 





is 


Streamlines 


en 

t 

a — 

Re =1.0*10° 


Vorticity Contours 


SSS 


a 97°70 
t 96.000 
rf 0.099 
Re =1.0*10° 





ei 


Streamlines 





78 


Streamlines 


Vorticity Contours 





qe, 


| Streamlines 


CX 
t OS. 
el 0.099 














rf = 0.099 
Re_=1.0*10" 






80 


streamlines 


Vorticity Contours 





ome 


streamlines 


Vorticity Contours 





82 


5.7 
t 4.560 
rf 0.150 
Re =1.0*10° 





83 


7a 
8.000 
0.150 
Re =1.0*10° 





84 





85 


a 
t 
rf 150, 
Re =1.0*10 





00 OQ1 0.2 0.3 0.4 0.5 0.6 0.7 08 0.9 10 





0.0 0.1 Ore Ons 0.4 0.5 0.6 0.7 08 0.9 10 
X/C 


86 





87 


24.3 

42.480} 
rf 0.150 
Re =1.0*10" 





88 


al 22.9 

t 45.920 
rf 0.150 
Re =10*10° 





PoeeOr “02 OdmmoO4s Ob 906 07 O08 09 1.0 
X/C 


89 


re 14.4 

t 57.520 
rf 0.150 
Re_=1.0*10" 





2.0 








MieeOeewOg 04 05 06 07 "08 a9 1.0 
GAC 


t | 
ri 0.150 
Re =1.0*10° 


sl 


0.0 


0.1 


OrZ 


0.3 


| 04 ¢ : 
X/C 


2) 


06 





io 


= 
— 


yUsloljyao) YT 


2 
°o 


- 
° 


1.2 


© 
ot 


>) ro) ro) 


yUusIDJa0g deiq 





0 


0.0 


0.1 


—~s-e wee ew ew a 


> > ee 
2° tae 


yUaToyjao) JUSWIOP, 





-0.5 


10.0 


5.0 


0.0 


a3 


Streamlines 


-Vorticity Contours 





94 


Streamlines 


Vorticity Contours 





2 


Streamlines 


oS aye. 
| 12.000 
ii GilaO 
Re =10*10' | 





0:6 


Streamlines 


Vorticity Contours 


X 


* 0.150 
me = 1.0*10 





oy. 


Streamlines 


15.870 
20.000 
Isl Oris 


i Vorticity. Contours 





98 


Streamlines 


Vorticity Contours 





a2 


Streamlines | 


OX 

t | 
raj 0.150 
Re =1.0*10° 





00 


Streamlines 


AOE 
32.000 


0.150, 
=1.0*10 





LOO 


Streamlines 


: _. se 
6 
Al 


OX 

t | 

rele 0.150 
Re =1.0*10° 





Or 


Streamlines 


ZZ E 
SS 


0 24 B49 
t 40.000 
Yr 


f 0.150. 
Re =1.0*10 


Vorticity Contours 


a 24.849 
t 40.000 
ei 0.150 
Re =1.0*10° 





Oe 


Streamlines 


i 


| Vorticity Contours _ 


aS 


a 23.748 
t 44.000 
mi = Oho 
Re =1.0*10° 





104 


paras 


Streamlines 


Vorticity Contours 





GS 


Streamlines 


a 18.911 
t 52.000 
et 0.150 
Re =1.0*10° 





106 


Streamlines 





O07 


Streamlines 





108 


streamlines 


OX 9.413 
t 64.000 
mire— iO 150 
Re =1.0*10° 





09 


i 
| | 
i 


or 
= 
SS 

© 


7017 ~ 
68.000 
0.150 


Streamlines 





=1.0+107 | © 


LI 


L 
ai 


Re =1.0*1 Q" | 


Streamlines 


i: Vorticity Contours 





JEIEME 


Streamlines 


Vorticity Contours 


on Olen 

76.000 
rf 0.150 
Re = 1.0*10° | 





PZ 


0.4 


X/C 


a3 


0.6 


a Sl 

t 1.040 
Jef 0.248 
Re = 1.0*10° 








ia 

4 880 
rf 0.248 
Re =1.0*10° 


93 

7.040 
rf = 0.248 
Re =1.0*10° 





0.0 0.1 0.2 0.3 0.4 0.5 06 On7 0.8 0.9 10 
X/C 


114 


el 

11.520 
fei 0.248 
Re =1.0*10° 





{Laks 





116 





LAT 


24.9 

21.920 
Ob = 024s 
Re =1.0*10° 


a = 244 

t = 25.520 
ri = 0248 
Re =1.0*10° 





Is: 


t | 
rf 0.248 
Re =1.0*10° 





0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
X/C 


ies, 


a Wee 

it 32.480 
ene 0245 
Re =10510" 





0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
X/C 


2.0 


X 

t | 

rf 0.248 
Re =1.0*10° 





0.0 0.1 0.2 0.3 0.4 0.5 0.6 OF 0.8 0.9 0 


2 





20.04 


15.04 


10.0 


—Cp 


5.0 





00 O1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 


2 


WUDIOIJJa0D WUT 





yUsIOIJa09g deig 


ow ow 


x 
= 


0 
, t 
t 
i] 
t 
f 8 
$ 
I 


° = ON 


SS 
-) 
| 


=073) 


yUSTOIJJaoy JUSUIOW 





f 
o 
] 


25.0 


0 


9.0 10.0 15.0 20 


0.0 


ne 


streamlines 


“Vorticity Contours 





124 


Streamlines 


Vorticity Contours 





leva: 


Streamlines | 


O 15.718 
tc ZOO 
rei O24 5 
Re =1.0*10° 





ZAG 


Streamlines 


Vorticity Contours 





le 7 


Streamlines 


= 
t 


a 24.195 
t 20.000 
rf 0.248 
Re =1.0*10° 





128 


OX 

t | 

ai 0.248 
Re =1.0*10° 


Streamlines 


Vorticity Contours 





IAS, 


Streamlines 


a an (OD 
t 28.000 
eli 0.248 
Re =1.0+10° 





red 








streamlines 







[ 


SSS SE 


fc = 18.260 
me 32-000 
im = 0.248 
Re =1.0*10° 


Vorticily Contours 


O 
t 





Jesal 


Streamlines 


t | 
eae 0.248 
Re =1.0*10" 





dies 2 


Streamlines 


a 8.089 
t 40.000 
esl OLeaus' 
Re =1,0*10° 


Vorticity Contours 





is 


Streamlines 





134 


AQAOXDANAAQAAAQAAADAAAANAAAANH 


APPENDIX B 


CODE LISTING 


PROGRAM GEOM 


TAPE1 oD 
TAPE2 00 
TAPES 00 
TAPES OD 
TAPE6 OD 
TAPE10 OD 


TAPE14 oD 
CALLS OD 


SSSSSHESSCSSSCSSVCSCesese Ss Ses se Sess Se sSesesesegecseseseseseseseeeteeaeseese 


SRHHHAHKS HERES KEKE REESE EERE KE KER HHENEHER EEE RENEE KES 
GEOMETRY PLOTTING PROGRAM — DISSPLA VERSION 
LAST REVISION 4—1-87 


PRINCIPAL INVESTIGATOR @D DR. J.C. WU 

AUTHORS @D MIKE PATTERSON, ISHMAEL TUNCER 
GEORGIA INSTITUTE OF TECHNOLOGY 
(404) 894-3028 


INPUT TO PLOT1 

input TO ZONST 

INPUT TO LOADS 
GENERAL INPUT 

GENERAL OUTPUT 

XYZ INPUT TO PLOT3D 
INPUT TO PLOT2, PLOTS 


NONE 


IMPLICIT REAL#8 (A-H,0-Z) 


PARAMETER 
PARAMETER 
PARAMETER 
PARAMETER 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENS ION 
DIMENSION 
DIMENSION 


(IDIM=80 , JDIM=60) 

(IP1=81, JP1=61) 

ice 
INOR1=20) 


UC(IP1,JP1),VC(IP1,JP1),HREAL(IDIM) ,HIMAG(JDIM) 
CS(KFC1, IDIM) ,SN(KFC1, IDIM) ,HSTAR(IDIM, JDIM) 
R1(JP1),R2(JP1),RiD(JP1),RP(JP1,KFC2),RL(JP1) 
H( IDIM, JDIM) , SGMA(KFC1) 

A2(JDIM) ,A4(JDIM) ,C2(JDIM) ,C4(JDIM) ,D5(JDIM) 
AS2(KFC1) ,BS2(KFC1) ,CS2(KFC1) ,DS2(KFC1) 
YN(KFC1,JDIM) ,CCP(KFC1,JDIM) 

IWK(JDIM, JOIM) , INOR(INOR1 , JDIM) ,COEF(IDIM, 2) 
DTT 101) omer JbIy ation TET 
XP(IP1,JP1),YP(IP1,JP1),XT(IDIM),YT(IDIM) 


COMPLEX W,W1,Z,2Z1,0ZDW, OWOT ,HIAMG ,HSTAR 


READ(5, #) 
ne ec) 
READ(5,* 


IM=IDIM 
JR=JDIM 
JR1=JP 1 
KFC=KFC1 


RE 
C1,CC2,DS 
CSQ,GAMA,SIGMA,AL 


PI=4.*ATAN(1.) 
DTET=2. #PI/FLOAT( IM) 


IM2=1M/2 


C..COMPUTE RADIAL GRID DISTRIBUTION IN COMPUTATIONAL PLANE 


DO 100 J=1,UR1 
S=FLOAT(J—1)#DS 
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REE (Sec 
R2(J )=EXP(S+C1+.5*DS)+CC2 
10@ CONTINUE 
DO 105 J=1,JR 
105 R1D(J)=R1 (U+1)-R1 (J) 


C..COMPUTE SCALE FACTOR OF TRANSFORMATION 


DO 116 J=1,JR 
DO 110 I=1,IM 

PHI=FLOAT(I-1) *DTET 

vpmecsavaceaian 

WI=R2(J)*SIN(PHI ) 

W=CMPLX (WR, WI) 

Z1=W+GAMA 

DZDW=1.-CSQ/(Zi*Z1) 

Hi=CABS(DZDW) 

H(1,J)=Hi eHi 

HREAL(I) = REAL(DZOW) 

HIMAG(J) = AIMAG(DZDW) 

HSTAR(I,J) = CMPLX(HREAL(1),-HIMAG(J)) 

110 CONTINUE 

WRITE(2@) HSTAR 


C..ARRAYS USED IN KINETICS 


DO 115 J=1,JR 
A2(J )=R1(J+1)/(DS*(R2( J )-CC2) ) 
A4(J )=-A2(J) eAL/(RE#DS*(R1(J+1)-CC2) ) 
C2(J )=R1(U)/(DS*(R2(JU )—CC2)) 
c4(J \=C2(J)*AL/(RE*DS*(R1(J)—CC2) ) 
D5(J J=AL/(R2(J) *RE*DTET*DTET) 
115 CONTINUE 


C..COMPUTE GEOMETRIC COEFFICIENTS FOR COMPUTATION OF VELOCITY 
C FOURIER COEFFICIENTS BY INTEGRAL RELATIONS 


DO 120 J=1,JR1 
aa Se 
RL(J)=LOG(R1(J)) 
DO 120 K=2,KFC+1 
RP(J,K)=RP(J,K-1) *RP(J,1) 
12@ CONTINUE 


C..COMPUTE COORDINATES AND DERIVATIVES OF SOLID SURFACE 
C FOR USE IN LOADS CALCULATIONS 


DO 125 I=1,IM 
PHI=FLOAT(I—1)*DTET 
CO=COS (PHI) 
SI=SIN(PHI ) 
W=CMPLX(CO,SI) 
DWDT=CMPLX(-SI,CO) 
Z1=W+GAMA 
Z=Z14+CSQ/Z1+SIGMA 
eee) 

YB(1 )=AIMAG(Z) 

DZDW=1.-CSQ/(Z1#Z1) 

ee OND) 

YT(1)=AIMAG(DZDW*DWDT ) 
125 CONTINUE 


C..COMPUTE CORRELATIONS OF VELOCITY COMPONENTS BETWEEN 
C INERTIA COORDINATE SYSTEM AND BODY COORDINATE SYSTEM 


PG 


DO 130 J=1,JR1 
DO 130 I=1,1M 
PHI=FLOAT(I-1) #OTET 
WR=R1(J)*COS(PHI) 
WI=R1(J)*SIN(PHI ) 
W1=WR+GAMA 
W2=W1 «W1+WI «WI 
XR=W14+CSQ«W1/W2+SIGMA 
YR=wI-CSQ=eWI /W2 
XC=1 .4CSQ* (WI «WI—-W1 #W1) /(W2*W2 ) 
XN=—CSQ¢2. «WI «W1/(W2«W2) 
WAR re ar vraxe 
VC(1,J)=YR*XN—XR*#XC 
13@ CONTINUE 
DO 135 J=1,JR1 
2 MRR oe ay 
VC(IM+1,J)=VC(1,J) 
135 CONTINUE 


C..COMPUTE COORDINATES OF VORTICITY GRID IN PHYSICAL PLANE 
C AND IN COMPUTATIONAL PLANE 


DO 140 J=1,JR1 
DO 140 I=1,IM 
PHI=FLOAT(I-1) #OTET 
Perot) <cineeut 
WI=R2(J)*SIN(PHI ) 
W=CMPLX (WR, WI) 
Z=W+GAMA+CSQ/ (W+GAMA )+S1GMA 
XP(1,J)=WR 
YP(1,J)=wI 
X(I ,J)=REAL(Z) 
Y(1,J)=AIMAG(Z) 
14@ CONTINUE 
DO 141 J=i,JR1 
ean eee eat 
YPC IM+1,J)=YP(1,J) 
141 CONTINUE 
WRITE(1@) IOIM,JDIM 
WRITE(10) | aay a RLS ae a ta 
1 (Y(I,J),1=1,IDIM) , J=1,JDIM) 
DO 145 I=1,IM 
PHI=FLOAT(I-1)*OTET 
CS(1,1)=1. 
SN(1,1)=@. 
CS(2,1)=COS(PHI ) 
SN(2,1)=SIN(PHI ) 
145 CONTINUE 
DO 150 K=3,KFC 
DO 150 I=1,IM 
L=1+MOD( (K-1)#( 1-1), 1M) 
CS(K,1)=CS(2,L) 
SN(K,1I)=SN(2,L) 
15@ CONTINUE 


C..COMPUTE FOURIER SERIES SMOOTHING FUNCTION 
DO 155 K=2,KFC 
A1=FLOAT(K—-1) #DTET 


SGMA(K)=SIN(A1)/A1 
155 CONTINUE 


Sy 


C..DETERMINE FOURIER COEFFS. OF BOUNDARY VORTICITY CONDITION 
C ON INERTIA COOR. (COEF(I,1) FOR RADIAL COMP., COEF(I,2) FOR 
C TANGENTIAL COMP. ) 


DO 160 K=1,KFC 
AS2(K)=0. 
BS2(K)=0. 
CS2(K)=0. 
0S2(K)=0. 
16@ CONTINUE 
DO 161 J=1,IM 
GOEL CUM egy Ue CS (2, Stee aay / FLOAT(IM2) 
COEF (I,2)=-(VC(1I,1)*CS(2,1)-UC(I,1)*#SN(2,1)) / FLOAT(IM2) 
161 CONTINUE 
DO 162 I=1,IM 
pele Saunier oD 
CS2(1)=CS2(1)+COEF(I,2) 
DO 162 K=2,KFC 
B82 (SAS eEnt  NSCSIOD 
BS2(K)=BS2(K)+COEF(I,1)*SN(K,1) 
082 (i face) COR Say 
DS2(K)=0S2(K)+COEF(I,2)*SN(K,I 
162 CONTINUE 
DO 163 K=2,KFC 
AS2(K)=AS2(K)#SGMA(K) 
BS2(K)=BS2(K) *SGMA(K) 
CS2(K)=CS2(K)*#SGMA(K) 
DS2(K)=DS2(K) #SGMA(K) 
163 CONTINUE 


C..COMPUTE ’YN’ 


DO 165 I=1,KFC 
DO 165 J=i,JR 
IS=] 
IF(I .GT. IM/3) THEN 
YN(I,J)#SQRT((X(1I,J)-XB(1))##2+(Y(1I,J)~YB(1))##2) 
ELSE 
DIF=x(1I,J)—-xB(1) 
IF(ABS(DIF) .GT. @.001) THEN 
INC=1 
IF(DIF .GT. @.) INC=1 
164 IS=IS+INC 
IF(IS .GT. KFC .OR. IS .LT. 1) THEN 
YN(I J )=Y(1,d) 
ELSE 
DIFN=X(I,J)-X8(IS) 
IF(DIFN*DIF .LT. @.) THEN 
ISM=IS—INC 
Y 1S=YB(1SM)+(X(I,J)-XB( ISM) )*(YB(1S)—YB( ISM) ) 
1 /(XB(1S)-XB( ISM) ) 
YN(I,J)=Y¥(I,J)-YIS 
SE 


GO TO 164 
ENDIF 
ENDIF 
ELSE 
YN(I, J)=Y(1,Ju)-YB(1) 
ENDIF 
ENDIF 
165 CONTINUE 


C..COMPUTE ’INOR’ 


IyBire: 


DO 170 Il=2,1M/4 

I=I] 

INOR(II,1)=I 

DO 170 J=2,UR 
INC=1 
IF(X(1,J) «LT. X€1I,1)) INC==1 
IF(ABS(X(II,1)—X(I+INC,J)) .LT. ABS(X(II,1)-X(I,J))) I=I+INC 
INOR(II,J)=1 

17 CONTINUE 


C..COMPUTE ° IWK’ 


DO 175 JN=1,JR-1 

I=INOR(2,JN) 

IWK (JN, JN)=I 

DO 175 J=JN+1,JR 
INC=1 
IF(Y(I,J) .GT. YCI,JN)) INC=1 
IF(ABS(Y(1,JN)-Y(I+INC,J)) .LT. ABS(Y(I,JN)-Y(I,J))) I=I+INC 
IWK(JN,J)=I 

175 CONTINUE 


C..CALCULATE GEOMETRIC COEFF. ON PRESSURE COMP. 


DO 180 J=1,JUR 
DO 179 I[=3,KFC 
179 CCP(I,J)= 1. / FLOAT(I-2)*(1./R1(J)*#(I-2) 
1 -1./R1(J+1)«e(I-2) ) 
oP (1. J} R1D(J) 
CCP(2,J)= LOG(R1(J+1)/R1(J) ) 
18@ CONTINUE 


C..INTEGRATE OVER AIRFOIL SURFACE TO COMPUTE AIRFOIL AREA; 
C ASSUMES AIRFOIL IS SYMMETRIC 


DPHI=P1/200. 

CO=—1. 

SI=0. 

W=CMPLX(CO,SI) 

Z1=W+GAMA 

Z=Z1+CSQ/Z1 

X1=REAL (Z) 

Y1=AIMAG(Z) 

PHI =P I 

Ai=0. 

DO 200 I=1,200 
PHI=PHI-—DPHI 
CO=COS( PHI) 
SI=SIN(PHI) 
W=CMPLX(CO,SI) 
Z1=W+GAMA 
Z=Z1+CSQ/Z1 
X2=REAL (Z) 
Y2=AIMAG(Z) 
A1=A1+.5*(Y2+Y1) #(X2-X1) 
X1=X2 
Y1=Y2 

200 CONTINUE 

AF=2. #A1 


C..OUTPUT TO THE VARIOUS FILES 


WRITE(1) AL,X,Y, INOR, IWK 
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WRITE(2) OTET,OS,AF,AL,RE 
WRITE(2) CS,SN,UC,VC,H 
WRITE(2) RP,RL,SGMA 
WRITE(2) A2,A4,C2,C4,05 
WRITE(2) AS2,BS2,CS2,0S2 
WRITE(2) R1,R2,R10,YN 
WRITE(2) INOR, IWK,CCP 


WRITE(3) AL,RE,OTET,XB,YB,XT,YT 
WRITE(3) CSQ,GAMA,SIGMA,C1,CC2 


WRITE(6, 14 

WRITE(6,4@) IM,JR,SIGMA,CSQ,C1,CC2,DS,GAMA,AL,RE 
WRITE(6,2@) (H(I,1),I=1,IM) 

WRITE(6,30) R1 


WRITE(14) CSQ,GAMA,SIGMA,AL 
WRITE(14) XB,YB,XP,YP 


1@ FORMAT(//,° INPUT TO GEOM’ ) 
20 FORMAT(//, ‘SQUARE OF THE SCALE FACTORS AT FIRST RING@D* 
1 ,//(1@E13.6)) 
3@ FORMAT(//,° GRID DISTRIBUTION IN R DIRECTION 6D '//(10E13.6)) 
4@ FORMAT(//,/1iX, 


1 ‘IM= ’,T10,15,/1X, 
2 ‘JR= ',71@,15,/1X, 
3 'SIGMA= °,T10,F10.7,/1X, 
4 °*CS@=.%,T10,F10.7./1x, 
5 ‘Ci= ',T10,F10.7,/1X, 

6 'CC2= °,T10,F18.7,/1xX, 
7 'DS= *,T10,F10.7,/1X, 
8 ‘GAMA= °,T10,F10.7,/1X, 
9 ‘AL= °,T10,F10.7,/1X, 

1 °RE= °,T10,F10.1) 

STOP 

END 
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18008000 


—4.927978847591 .99275887608358 .12935603008237 


.8061367 -.0647384 .9170954997604 3.619058974007 
INPUTS: NACA 9212 REVISED PARAMETERS 

RE 

C1,C2,D0S 


CSQ,GAMA, SIGMA, AL 
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PROGRAM ZONST 


PTET ETICTOCCTOCSCTOSCOCSOSCTOCOCOSCTOCSSOCC COC COT TEST eee eee Teese er eee 
2—D INCOMPRESSIBLE NAVIER-STOKES SOLVER — DISSPLA VERSION 
LAST REVISION@D 5-20-87 


PRINCIPAL INVESTIGATOR @D DR. J.C. WU 

AUTHORS @D MIKE PATTERSON, ISHMAEL TUNCER 
GEORGIA INSTITUTE OF TECHNOLOGY 
(404) 894-3028 


TAPE2 OD OUTPUT FROM GEOM 

TAPE4 00 OUTPUT FOR LOADS 

TAPES O00 GENERAL INPUT 

TAPE6 6D GENERAL OUTPUT 

TAPE7 OD INPUT FROM PREVIOUS RUN 
TAPE8 OD OUTPUT FOR NEXT RUN 
TAPES OD OUTPUT FOR PLOT2 

TAPE11 OD Q INPUT FOR PLOT3D 


CALLS 9D KMTCS 
KNTCS 
CPVAL 


SCESCSCSSESCKESCSCSSPSECK SSS SHSSE See SCSHSSSeeese Sse Sees esses esheeeeesese 


PARAMETER (IDIM=80, JDIM@60) 
PARAMETER (IPi=81,JP1=61) 
PARAMETER (KFCi=41,KFC2=42) 
PARAMETER (IDUM=4295, INOR1=20) 


DIMENSION AS2(KFC1),8S2(KFC1) ,CS2(KFC1) ,DS2(KFC1) 
COMMON/IO/LOOP , NT, NTMAX , NTOUT 

COMMON /DKIN/A2(JDIM) ,A4(JDIM) ,C2(JDIM) ,C4(JDIM) ,D5(JDIM) 
COMMON /VLB/VORLB , NLB 

COMMON/RGRD/R1D(JP1),R2(JP1) 

COMMON/SMT /SGMA (KFC1 ) 
COMMON/WFC/AA(JDIM,KFC1),BB(JDIM,KFC1) 

COMMON/RHS /WOW(IP1,JP1),POP(IP1,JUP1),WP(IP1),WB(IP1),R1(JP1), 
1 DUMM( IDUM) 

COMMON/ABCDS/AS1(KFC1) ,BS1(KFC1),CS1(KFC1),0S1(KFC1) 
COMMON/COR/RP(JP1,KFC2),RL(JP1) 

COMMON/DELTA/DS , DTET , DT 

COMMON/UNF /URR , URF , DFMX , NCC 

COMMON /VTEXT/KIN(IDIM) ,KEN( IDIM) 
COMMON/TRIG/CS(KFC1, IDIM) , SN(KFC1, IDIM) 
COMMON/SCALE/H(IDIM, JDIM) 
COMMON/VV/VOR( IDIM, JDIM) , VOROLD( IDIM, JDIM) 
COMMON/VEL/U(IDIM, JP1),V(IDIM, JP1) ,HSTAR(IDIM, JDIM) 
COMMON/TUR1/USTAR(IDIM) , EDDY(IDIM,JDIM) , IOTUR(IDIM), 

1 YN(KFC1,JDIM), IOT, IOTB 
COMMON/TUR2/INOR( INOR1,JDIM) , IWK(JDIM, JDIM) 
COMMON/VEH/UC(IP1,JP1),VC(IP1,JP1) 
COMMON/ZNS/IB1,1V1,1V2,1B2,IV1R 

COMMON/GRD/IM, IM2,KFC,JR,N 

COMMON/COE/AF ,UI,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 
COMMON/CPC /CCP(KFC1, JDIM) ,CP(IP1) 

COMMON Q1(IDIM, JDIM) ,Q2(IDIM, JDIM) ,Q3(IDIM, JDIM) ,Q4(IDIM, JOIM) 
COMPLEX HSTAR 


REAL ALLP(2) 
DATA ALLP/ 19.5,20.0/ 
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REWIND 
REWIND 
REWIND 
REWIND 
REWINO 
REWIND 
REWINO 
REWIND 20 


MNQAQAQAQNQN 
OOnN AON 


C..READ GENERAL INPUTS AND ECHO THEM 


READ(5,+*) ICST,RF,ALPS, ICTUR 
READ(5,*) WMIN,OFMX,DRMX,KMAX,NCC 
READ(5,*) URBI,URBP,URR 

READ(5.*«) IV1,1B1,1B2,1V2,NPL,NLB 
ee DTI ,DTINC, NTMAX 
READ(S,<) NTPL,NTOUT,NTLO 


WRITE(6,54) 

WRITE(6,*) ICST,RF,ALPS, ICTUR 
WRITE(6,*) WMIN,OFMX, ORMX , KMAX , NCC 
WRITE(6,*) URBI,URBP,URR 
WRITE(6,*) IV1,1B1,1B2, 1V2,NPL,NLB 
WRITE(6,*) OTI,DTINC,NTMAX 
WRITE(6,*) NTPL,NTOUT,NTLO 


C..INPUTS FROM GEOM 


READ(2) OTET,DS,AF,AL,RE 
READ(2) CS,SN,UC,VC,H 
READ(2) RP.RL,SGMA 
READ(2) A2,A4,C2,C4,D5 
READ(2) AS2,BS2,CS2,DS2 
READ(2) R1,R2,R1D,YN 
READ(2) INOR, IWK,CCP 
READ(2@) HSTAR 


NALLP=1 

IM=IDIM 

JR=JOIM 

KFC=KFC1 
URF=1.—-URR 
NLBP=NLB+1 

P[=4. sATAN(1.) 
IMi=IM+1 
1ViR=IM+IV1 
IM2=IM/2 
RDTET=DTET *R1 (NLBP) 
VSC=AL/RE 

FR=2. #RF/AL 
ALP=ALPS*PI/180. 
OMG=2 . 


OMGD=0 . 
C..START INITIAL SOLUTION OR READ PREVIOUS ITERATION RESULTS 


IF(ICST.LT.2) THEN 
DO 100 K=1,KFC 
AS1(K)=0. 
BS1(K)=0. 
CS1(K)=0. 
DS1(K)=0. 
108 CONTINUE 
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ENDIF 
IF(ICST.EQ.@) THEN 
NT=0 
IF(OTI.EQ.0.) DT=0.08 
T=0.0 
N=15 
VORLB=@ . @ 
UI=COS(ALP) 
VI=SIN(ALP) 
DO 165 I=1,IM 
KIN(I)=N 
U(1,1)=@. 
V(J,1)=@. 
DO 105 J=1,JR 
VOR(I ,J)=0. 
105 CONTINUE 
AA(1,1)=@. 
AA(1,2)=2.*VI/(RP(2,1)—RP(1,1)) 
BB(1,2)=-2.*U1/(RP(2,1)—-RP(1,1)) 
C..POTENTIAL FLOW SOLUTION 


CALL KMTCS 
DO 106 I=1,IM 
106 VOR(I, 1)=(AA(1,2) *CS(2, 1)4+8B(1,2)*SN(2,1))/H(I,1) 
ENDIF 


C..READ PREVIOUS RUN AND RESET TIME IF A NEW MOTION IS SPECIFIED 


I=ICST 
IF(ICST.GE.1) READ(7) ICST,NT,N,.KIN,T,DT,VOR,U,V,VORLB 


IF(1.GT.1CST) THEN 
IF(I.NE.1) THEN 
WRITE(6,51) 

NT=@ 
T=@.Q 
ENDIF 

ENDIF 

ICST=] 

IF(DTI.NE.@.) DT=0TI 


COSS/ILITI/ TIME STEP LOOPS///////IIII/L/111// 
C..START COMPUTATIONS FOR SUBSEQUENT TIME STEPS 
C FROM TIME LEVEL NS1 TO NTMAX (DO LOOP 10@1) 


DO 1001 LOOP=1,NTMAX 
TO=T 
NTO=NT 
DTO=DT 
DT=DT+DTINC 
IF (DT.GT..08) DT=.e8 
T=T+DT 
NT=NT+1 
DO 11@ J=1,N 
DO 11@ I=1,IM 
11@ VOROLD(I,J)=VOR(I,J) 


C..DEFINE AIRFOIL MOTION 
IF(ICST.£Q.2) THEN 
ALP=(T-SIN(7.5*T)/7.5)*#FR#O.5 


OMG=—(1.-COS(7.5#T))*FReO.5 
OMGD=—SIN(7.5*#T)*FR*@.5*7.5 
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ENDIF 
IF(ICST.EQ.3) THEN 
ALP=FReT 
OMG=—FR 
IF(ALP.GT.@.6) THEN 
ALP=0.6 
OMG=0. 
ENDIF 
ENDIF 
IFC ICST.EQ.4) THEN 
ALP=(15.-10. sCOS(FReT))*#P1/180. 
OMG=-10.*FReSIN(FR*T)*P1/180. 
OMGD=—10. #FR*FR*COS(FR*T) sP1/180. 
ENDIF 


C..VELOCITY BOUNDARY CONDITION ON BODY 


IF(ICST.GT.1) THEN 
DO 115 K=i,KFC 
AS1(K)=AS2(K) #OMG 
BS1(K)=BS2(K ) OMG 
CS1(K)=CS2(K) OMG 
DS1(K)=DS2(K) #OMG 
115 CONTINUE 
ENDIF 
ALPD=ALP #180. /PI 
UI=COS(ALP) 
VI=SIN(ALP) 
IF (INT(ALPD*10.) .€Q. INT(ALLP(NALLP)*10.)) THEN 
ICPL =@ 
ICLD = @ 
ICOUT = @ 
NALLP = NALLP+1 
WRITE(6,61) NT,T,DT,ALPD,OMG,OMGD 
ELSE 


C PONT NTEL) 
C ICLD=MOD(NT,NTLO 

C ICOUT=MOD (NT ,NTOUT) 

C WRITE(6,61) NT,T,OT,ALPD,OMG,OMGD 
c 


.. SOLVER FOR KINETICS 


CALL KNTCS(WMIN,URBI ,URBP ,KMAX, DRMX , KODE, KKK, ALP) 
IF(KODE.EQ.1) GO TO 5800 


C..UPDATE ARRAY KIN 


DO 120 I=1,IM 
JI=KEN(1)+3 
IF(JI.GT.JR) JI=JR 
IF(N.LT.JI) NeJI 
KIN(I)=vJI 

12@ CONTINUE 

IF(N.LT.NPL) N=NPL 


C..SOLVER FOR KINEMATICS 


CALL KMTCS 
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C..EVALUATE VORTICITY LEAVING BOUNDARY 


DO 125 I=1,IM 
IF(KEN(1).GE.NLB) VORLB=VORLB+DT *RDTET*(V(I ,NLBP) *@.5*(VOR(I,NLB) 
1 +VOR(I,NLBP) )-VSC*#(VOR(I,NLBP)—VOR(I,NLB) )/(R2(NLBP)—-R2(NLB) ) ) 
125 CONTINUE 


C..EVALUTE THE TOTAL VORTICITY INCLUDING SOLIO ROTATION 


WSUM1=2 . #OMG*AF+VORLB 
DO 13@ J=1,NLB 
ACA=R2(J) #DTET*R1D(J) 
DO 13@ I=1,IM 
WSUM1=WSUM1+VOR(I,J)*H(I,J) eACA 
13@ CONTINUE 
WRITE(6,57) KKK,WSUM1,VORLB 


C..CALCULATE CP 


IF(ICLD.EQ.@ .OR. ICPL.£Q.@ .OR. LOOP.£Q.NTMAX) THEN 
CALL CPVAL 
DO 135 I=1,1M 
IF(ICTUR.EQ.1 .AND. JOTUR(I).£Q.1) THEN 
WB(1I)=USTAR(I) #USTAR(1)/VSC 
IF(VOR(I,1).LT.@.) WB(1)=-WB(1) 


LSE 
WB(I)=VOR(I,1) 
ENDIF 
135 CONTINUE 
WB(IM1 )=WB(1) 
WRITE(4) NT,NTPL,T,RF,WB,CP,OMG,OMGD,ALP 
ENDIF 


C..OQUTPUT AT EVERY NTOUT TIME STEP 
IF(ICOUT.EQ.0 .OR. LOOP.EQ.NTMAX) THEN 


c Mle 8765) (vor (Fae oe 
c WRITE(6, 67) a 
é WRITE(6,69) (KEN(I),I=1,1M) 
c WRITE(21) W1,W2 
C IF(ICTUR.EQ.1) WRITE(6,70) IOT,10TB,IM 
ENDIF 
C..CALCULATE STREAM FUNCTION AND VORTICITY AT VELOCITY 
C GRID POINTS AND STORE ON TAPE FOR GENERATING PLOTS. 
C INTEGRATION IS FIRST ORDER TRAPAZOIDAL RULE. 


IF(ICPL.EQ.0) THEN 
DO 14@ I=1,1M1 
140 POP(!,1)=0. 


C..INTEGRATE TANGENTIAL VELOCITIES 


DO 15@ I=1,I1M 
DO 149 J=2,NPL 
wow(1,J)=.5*#(VOR(1,J)+VOR(I,J-1)) 
POP(I,J)=POP(I ,J—1)-—.5*(U(CI,J—1)4+U(1,J))*R1D(J—-1) 
149 CONTI NUE 
wowW(1I,1)=2.*VOR(I, 1)—WOW(1,2) 
15@ CONTINUE 
DO 151 J=1,NPL 
WOW(1M1,J)=WOW(1,J) 
POP(IM1,J)=POP(1,J) 
151 CONTINUE 
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WRITE(9) NT,NPL,T,ALPO,RF,RE,wWOW,POP 
ENDIF 


1001 CONTINUE 
CSIIIIITIL11//////ENO OF TIME STEPSSS///////L///1111//11/1 


C..ABNORMAL EXIT. WRITE DATA FROM LAST COMPLETED TIME STEP 


5@0@ IF (KODE.£Q.1) THEN 

WRITE(6,52) KKK 
WRITE(6,53) NTO 
WRITE(6,63) (VOR(I,1),I=1, 1M) 
rie (6, 87) (u(1,2), 1=1, 1M) 
WRITE(6,69) (KEN(I),[=1, IM) 
IF(ICTUR.EQ.1) WRITE(6,70) IOT,IOTB,IM 
WRITE(8) ICST,NTO,N,KIN, TO,DTO, VOROLD,U,V,VORLB 
STOP 

ENDIF 


C..NORMAL EXIT. WRITE THE REQUIRED DATA FOR NEXT RUN 


WRITE(6,53) NT 
WRITE(8) ICST,NT.N,KIN,T,0T, VOR,U,V,VORLB 
WRITE(11) IDIM,JDIM 
WRITE(11) 1.,ALPD,1000000. ,1 
DO 600 I=1,IDIM 
DO 602 J=1,JDIM 
Qi(I,J)=1 
600 CONTINUE 


WRITE(11) ((Q1(I,J),1=1, IDIM) ,J#1,JDIM), 
1 Q2(1,J),1=1, IDIM) ,J=1,JDIM), 
2 Q3(1,J),1=1,IDIM) ,J=1,JDIM), 
3 Q4(1,J),1=1, IDIM) , J=1,JDIM) 

STOP 


51 FORMAT(//20X,’*AIRFOIL MOTION HAS BEEN CHANGED, TIME RESET*',//) 
52 FORMAT(/10X,°sNO CONVER. SOLU. WITHIN’, 15, 3X,’ ITERATIONS®’ Ee 
53 FORMAT(//,30X,’$$ NEXT RUN WITH NSTART =',16,° $$°) 

54 FORMAT(//3X,° INPUT FILE TO ZONST’'/) 


57 FORMAT(12X,ITKNT =',14,9X,’TVOR = °,£13.6,° VORLB =’,£13.6) 
61 FORMAT(/2X,9@(1H-)/1X," NT =",14,',',’ T =',F8.4,' DT =",F8.6 
AA =',F7.4, OMG =’,F6.4,” OMGD =',F6.4/2X,90(1H-)) 


63 FORMAT(//1H ,'VORTICITY ——FIRST RING——' 
1 ,'CCW FROM TRAILING EDGE@D’/, /(10E13.6)) 
67 FORMAT(/1H ,° TANGENTIAL VELOCITY ——J=2——' 
1 ,’CCW FROM TRAILING EDGE@D’ /, /(10E13.6)) 


69 FORMAT(/1H ,’VORTICITY EXTENT oD eee KEN ’ 
1 ,’ee08’ ,/,/(1H ,4(2X,1013))) 

70 FORMAT(/1H ,' TURBULENT REGION OD °,13,’- 1 AND °‘,13,'—',13) 
END 
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DMDOANQIDIAIAAAYNANO 


MOOAdAY 


SUBROUTINE KNTCS(WMIN,URBI , URBP , KMAX , DRMX , KODE, KKK ,ALP) 


SEECeeeeeeSSKKKRKSEKRKESRERSEESSEKCKRSEKCEKEKEESEKEEKEEKEESEEREKEKREREKESEE 


KINETICS OF THE PROBLEM 


CALLS OD WSURFT  WSURF 
EDDYS FOCFT 
VORTY 


CALLED BY @D ZONST 


SSSSSSKESSSEESRESSEESEKSESSESEHSESSESKEESKSEHESSESEKEKSSEKESKKSETESSEKRE SE HSS 


PARAMETER (1D IM=80, JDIM=60) 
PARAMETER (IP1=81,JUP1=61) 
PARAMETER (KFC1=41) 


COMMON/10/LOOP , NT , NTMAX , NTOUT 

COMMON /VTEXT/KIN(IDIM) ,KEN(IDIM) 

CORON S247 GAL IDIM JOE) ,GB(IDIM, JDIM) ,GC(IDIM, JDIM) 
COMMON/RHS/GD(IDIM, JOIM) ,DIP1(IDIM, JDIM) ,DIM1(IDIM, JDIM) 
COMMON /WFC/AA(JDIM,KFC1) ,BB(JDIM,KFC1) 

COMMON/COE/AF ,UI,VI,OMG,VSC,NPL,ICTUR, ICST, ICPL 
COMMON/TRIG/CS(KFC1, IDIM) , SN(KFC1, IDIM) 

COMMON/DELTA/DS , OTET, DT 

COMMON/ZNS/IB1,1V1,IV2,1B2,IViR 

COMMON/TUR 1/USTAR(IDIM),EODY(IDIM, JDIM) , IOTUR(IDIM), 

1 YN(KFC1,JDIM) , OT, IOTB 

COMMON/GRD/IM, IM2,KFC,JR,N 

COMMON/SCALE/H( IDIM, JDIM) 
COMMON/VV/VOR(IDIM, JDIM) , VOROLD(IDIM, JDIM) 
COMMON/VEL/U(IDIM, JP1), V(IDIM, JP1) 

CONG mene ABS Sal ta 

COMMON/DKIN/A2(JDIM) ,A4(JDIM),C2(JDIM) ,C4(JDIM) ,D5(JDIM) 
COMMON /FM/GAM(IDIM, JDIM) 


KODE=@ 

IKS=1 

1 CALE Eee) IKS=IM2+2 

IF(MOD(NT,NTOUT).£Q.@ .OR. LOOP.EQ.NTMAX) WRITE(6, 10) 


. COMPUTE FRICTION VELOCITY AND EDDY VISCOSITY 


IF(ICTUR.EQ.1) THEN 
CALL WSURFT 
CALL EDDYS 

ENDIF 


. CONSTRUCT DATA ARRAYS THAT ARE NEEDED IN COMPUTATION OF 
FINITE OIFFERENCE FORM OF VORTICITY TRANSPORT EQUATION 


TIME TERM —- FORWARD DIFFERENCE 
DIFFUSION TERMS = CENTRAL DIFFERENCES 
CONVECTION TERMS ~— 2ND UPWIND DIFFERENCES 


DO 101 I=IKS,IM 

IPP1=141 

IM1=I-1 

eae IPP1=IPP1—IM 
IF(IM1.LT.1) IM1=IM+IM1 
JL=KIN(I) 


C..DETERMINE IF REGION IS BL OR NS 


ICBL=0 
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IF((1.GT.IV1.AND.1.LT.1B1).OR.(1.GT.1B2.AND.1.LT.1V2)) ICBL=1 
C..CONSTRUCT MATRIX COEFFICIENTS FOR INTERIOR VORTICITY SOLUTION 


DO 100 J=2,JL 
TI=H(I,J)*#R2(J)/OT 
GA(1,J)=A4(J) 
GB(1,J)=T1-A4(J )-C4(J) 
GC(1I,J)=C4(J) 
GD(1,J)=T1i*VOROLD(I1,J) 


C..TURBULENT FLOW TERMS AND BOUNDARY LAYER TERMS 


IF(ICTUR.EQ.1 .AND. IOTUR(I).EQ.1) THEN 

JPPi=J+1 

IF(JPP1.GT.JR) JPP1=JR 

GA(I,J)=GA(1,J)+A4(J)*EDDY(1, JPP1) /VSC 

Se le —(A4(J)+C4(J)) *EDDY(1,J)/VSC 

GC(I,J)=GC(1,U)+C4(J) *EDDY(1, J-1)/VSC 

ENDIF 
IF(ICBL.£Q.1) THEN 
te a 
DIM1(I,J)=0. 
ELSE 

PA coed 

DIM1(I,J)=05(J 

GB(1,J)=GB(1,J)+2. «05(J) 

IF(ICTUR.EQ.1 .AND. IOTUR(I).EQ.1) THEN 
ee eare IE v)e08(J)-EDOYS IPPs 4) /vsC 
DIMi(I,J)=DIM1(I,u)+D5(J)*EDDY(IM1,J)/VSC 
GB(1,J)=GB(1,J)+2. *05(J) *EDDY(I,u)/VSC 

ENDIF 

ENDIF 


C..DISCRETIZATION OF CONVECTION TERM IN RADIAL DIRECTION 


VReV( I, J+1) 

VL=V(I,J) 

A2VR=A2(J)*VR 

C2VL=C2(J)*VL 

IF(VR.LT.@.) THEN 
GA(1,J)=GA(I1,J)+A2VR 
LSE 


GB(1, J)=GB(I,J)+A2VR 
ENDIF 
IF(VL.LT.@.) THEN 
GB(1,J)=GB(I,J)+C2VL 
E 


ELS 
GC(1,J)=GC(1,J)+C2Vt 
ENDIF 


C..DISCRETIZATION OF CONVECTION TERM IN TANGENTIAL DIRECTION 


UL=0.25*(U(IM1, J+1)+U(IM1 ,J)+U(I ,J+1)4+U(1,J)) 
UR=0.25«(U(I ,J+1)+U(1,J)+UCIPP1,J+1)+UCIPP1,J)) 
OTUL=UL/DTET 
OTUR=UR/DTET 
IF(UL.LT.@.) THEN 
GB(1,J)=GB(I1I,J)-DTUL 
ELSE 
DIM1(1,J)=DIM1(I,J)4+DTUL 
ENDIF 
IF(UR.LT.@.) THEN 
DIP1(1I,J)=DIP1(1I,J)-DTUR 
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EESE 
GB(1,J)=GB(1,J)+DTUR 
ENDIF 
188 CONTINUE 


C..NEUMANN TYPE B.C 


IF(JL.GE.JR) GB(I,JL)=GB(I,JL)+GA(I, JL) 
101 CONTINUE 

KKK=@ 

[CMR=2 

DMAXP=100. 

URB=URBI 


C..VORTICITY CONVERGENCE LOOP 


50@ CONTINUE 
KKK=KKK+1 


C..BOUNDARY CONCITION FOR KINETICS IN TURBULENCE 


IF(ICTUR.EQ.2) THEN 
DO 11@ I=1,IM 
IF(IOTUR(1).£Q.1) THEN 
IN=I 


IF(IN.GT.KFC) IN=IM-1+2 

YP2=YN(IN,2)*#USTAR(I)/VSC 

IF(YP2.GT.5.) THEN 
JOTUR(I)=@ 

ELSE 
VOR(I,2)=VOR(I,1) 

ENDIF 

ENDIF 
11@ CONTINUE 
ENDIF 


C..SOLVE FOR EACH ZONE; USE ONLY HALF OF GRID IF SYMMETRIC FLOW 


IF(ALP.NE.@.) THEN 
CALL VORTY(1B1,1B2,1,1ICS,@,KODE) 
IF((IV1+1).LT.1B1) CALL VORTY(IB1-1,IV1,-1,1CB1,1,KODE) 
IF(IV2.GT.IB2) CALL VORTY(1B2+1,1V2,1,ICB2,1,KODE) 
CALL VORTY(IV2,IViR,1, ICV,@,KODE) 

ELSE 
CALL VORTY(IM2+2,1B2,1,1CS,@,KODE) 
CALL VORTY(IB2+i, 1V2,1,1CB2,1,KODE) 
CALL VORTY(IV2,1M,1,ICV,@,KODE) 
DO 12@ I=2, 1M2 
1I=IM+2—] 
DO 120 J=2,N 

12@ VOR(I,J)=-VOR(II,J) 
ENDIF 
IF(KODE.GT.@) RETURN 


C..DETERMINE VOR. EXTENT AT EACH RADIAL LINE IN N-S REGION 
C AND STORE IN ARRAY KEN 


DO 130 I=1,IM 
DO 129 J=N,2,-1 
129 IF(ABS(VOR(I,J)).GT.WMIN) GO TO 130 
13@ KEN(I)=J 


C..DETERMINE FOURIER COEFFICIENTS OF VORTICITY AND EVALUATE 
C SURFACE VORTICITY BY UNDER—RELAXATION TECHNIQUE 


S60 


DO 140 I[=1,I1M 
JL=KEN(1) 
JLP=JL+1 
DO 139 J=1,JL 
139 GAM(I,J)=H(1,J)*VOR(I,J)/FLOAT(1IM2) 
IF(JLP.LE.N) THEN 
DO 138 J=JLP,N 
138 GAM(I,J)=0. 
ENDIF 
14@ CONTINUE 
CALL FOCFT 
CALL WSURF 
MI=@ 
MW I= 
WMAX=0 . 
DMAX=9 . 
Wi=0. 
DO 150 I=1,IM 
WizAA(1,1)*.5 
DO 149 K=2,I1MZ2 
149 Wi=W1+AA(1,K)*CS(K,1)+BB(1,K)*SN(K, 1) 
W1=W1+AA(1,KFC)*CS(KFC,1)*.5 
Wi=W1/H(1,1) 
URI=URB*H(I,1)**URBP 
WWeW 1 sUR1+(1.-UR1)*VOR(1,1) 
AWW=ABS (WW) 
IF(AWW.LT.10.) AWW=12. 
DW=ABS (WW-VOR(I,1))/AWW 
IF(DW.GE.DMAX) DMAX=DW 
IF (DW.GE.DMAX) MI=I 
IF (AWW.GE.WMAX) WMAX=AWW 
IF (AWW.GE.WMAX) MWI=] 
VOR( I, 1)=WW 
15@ CONTINUE 


C..ADJUST UNDER-RELAXATION PARAMETER URB. 
C EXIT KNTCS IF CONVERGED, CONTINUE ITERATIONS IF NOT. 
C ABORT IF MAXIMUM ITERATIONS EXCEEDED. 


IF(KKK.EQ. ICMR) THEN 
DDM= (DMAXP—DMAX ) /DMAX 
IF(DDM.LT.@.) URB=0.90*URB 
IF(DDM.LT.@.1.AND.DDM.GT.@.) URB=1.@8sURB 
IF(DDM.GT.@.1) URB=1.05*URB 
IF(URB.GT.@.65) URB=@.65 
ICMR=ICMR+2 

ENDIF 

DMAXP=DMAX 

IF(MOD(NT,NTOUT).EQ.@ .OR. LOOP.EQ.NTMAX) 

+ WRITE(6,12)KKK,DMAX,MI,VOR(MI,1),WMAX,MWI,URB, ICS, 1CB1, 1CB2, ICV 
IF (DMAX.LE.DRMX) RETURN 
IF(KKK.LT.KMAX) GO TO 502 


C..NO CONVERGENCE OCCURED; RETURN TO MAIN WITH KODE=1 


KODE=1 
IF(MOD(NT,NTOUT).NE.@ .OR. LOOP.NE.NTMAX) THEN 
eue (S12) 
WRITE(6,12)KKK,DMAX,MI , VOR(MI,1),WMAX,MWI,URB, ICS, 1CB1, 1CB2, ICV 
ENDIF 
RETURN 
1@ FORMAT(//,25X, SURFACE VORTICITY INFORMATIONOD’ , 20x, 
1 ‘INTERIOR VORTICITY INFORMATION@D’/’ ITER’ ,6X, 


lige ge 


2 ’REL. MAX. I VOR(I,1)’,9X, MAX. VAL. I URB 

3 9X,’ITSED ST,BL1,BL2,NS’/,13X, 'DIFF.’,29X,’OF VOR.’) 

12 FORMAT(14, 4%, £135.6, 14, 1X,E135.5,6X%,611.5, 14, 69.2, 16x elo) 
END 


La 


> 


SUBROUTINE WSURFT 


SSSSESKRSH SSS SSSSKKS SSK SE KSSESSAHEARHESSSESESSSESSKEESKSEKREKSKEKEKEKRRESESE 


EVALUATE BOUNDARY VORTICITIES IN TURBULENT FLOW REGION 
CALLS @D0 USTARF 
CALLED BY 0D KNTCS 


SSKSASESCSS HSS SSSR KSKCAHKKSSKSKASHKASSKS SKS SRSKSHRSEKRKEKARSSKSKKESK ASKER SBE 


AAAAAAAANQN 


PARAMETER (IDIM=80, JD IM=60) 
PARAMETER (KFC1=41) 


COMMON/SCALE/H( IDIM,JDIM) 
COMMON/VV/VOR(IDIM, JDIM) , VOROLD( IDIM, JDIM) 
COMMON/TURI/USTAR(IDIM) ,EDDY(IDIM,JDIM) ,IOTUR(IDIM), 
1 YN(KFC1,JDIM), 1OT, IOTB 

COMMON/COE/AF ,UI,VI,OMG,VSC,NPL,ICTUR, ICST, ICPL 
COMMON/GRD/IM, IM2,KFC,JUR,N 


C..EVALUATE USTAR 


OMG2=2. *OMG 

USTAR(1)=0. 

DO 100 I=2,IM 

IN=] 

IF(IN.GT.KFC) IN=IM-I+2 

GM=(VOR(I,1)—OMG2)*YN(IN, 1) 

AGM=ABS (GM) 

YP=SQRT ( AGM* YN(IN,1)/VSC) 

USTAR(I )=AGM/YP 

IF(YP.GT.5.) USTAR(I)=USTARF(USTAR(I),AGM, YN(IN,1),5.,-3.@5, YP) 

IF(YP.GT.3@.) USTAR(I )=USTARF(USTAR(I),AGM,YN(IN,1),2.5,5.5, YP) 
10@ CONTINUE 

RETURN 


END 
FUNCTION USTARF(USTR,UV,YN,A,C.YP) 


SVSSKSCSCSKS SSCS SH SHS SHS SS SKC KCC SCS SP SCC SCS ses swC SCS ss se eseseseeseeseseseeeaeeae gs 


EVALUATE USTAR ITERATIVELY IN INERTIAL LAYER 
CALLS OD NONE 
CALLED BY @D WSURFT 


SIRSKSKC SHAH SSSKSK SH SSKSASSKSSSEKHKKCSHTSESSesesee#wesese see seseseseseeaeggegeseegee 


AAIANAAAAAAQA 


COMMON/COE/AF ,UI,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 
DO 10@ I=1,20 
YP=YN#USTR/VSC 
USTARF=UV/(A*ALOG(YP)4C) 
ABSER=ABS( (USTARF—USTR) /USTARF ) 
IF(ABSER.LT.0@.005) RETURN 
USTR=USTARF 
18 CONTINUE 
WRITE(6,1@) YP, USTR,ABSER 
1@ FORMAT(2X,°NO CONVERGENCE IN USTARF. YP,USTAR,ABSER@D’ ,3F9.5) 
STOP 
END 
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MOOQOOdNdOQNAAIAIO 


SUBROUTINE EDDYS 


SRHESSSRREERSSSHESSHSSEKSSEHSSSERAEHESHESSSEHESSRESSESSERHERSERSESERSEASE KRESS 


EVALUATE EDDY VISCOSITY 
CALLS OD NONE 
CALLED BY @D KNTCS 


SESKSEEPEESHHKSTCRKEKESCHMERSSEHRHEKCATSHESSHSHRSSHS SSS KSKESseeSeeseeesveesee s 


PARAMETER (IDIM=86, JDIM=60) 
PARAMETER eee aD 
PARAMETER (KFC1=41) 
PARAMETER (INORi1=20) 


COMMON/SCALE/H( IDIM, JDIM) 

COMMON/10/LOOP ,NT , NTMAX , NTOUT 

COMMON/VV/VOR(IDIM, JDIM) , VOROLD(IDIM, JDIM) 

COMMON/TUR1/USTAR(IDIM) , EDDY (IDIM, JDIM) , IOTUR(IDIM), 
YN(KFC1,JDIM), 10T, IOTB 

COMMON/TUR2/INOR(INOR1, JDIM) , IWK(JDIM, JDIM) 

COMMON/GRD/IM, IM2,KFC,JR,N 

COMMON/COE/AF ,UI,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 

COMMON/VTEXT/KIN(IDIM) ,KEN(IDIM) 

COMMON/VEL/U(IDIM, JP1),V(IDIM, JP1) 


IMB=IM+2 
DO 185 IJ=2,1M 
J=[] 
IOTUR(I )=6 
PMAX=@. 
YMAX=@ . 
VDMX=@ . 
JL=KING I) +1 
IF(JL.GT.JUR) JL=JR 
IFGetCLs3} 67-3) JL=KIN(1+1) 


IF(KIN(I-1).GT.gL) JL=KIN(I-1) 


C..CALCULATE TERMS INDEPENDENT OF Y AND INNER VISC. 


C 


INCLUDE CORRECTION FOR THE NORMAL DIRECTION 


DO 10@ J=2,JL 
IF(II.LT. INOR1+1 ) ce ies 
IFCII.GT.JP1) I=sIMB—INOR(IMB=II,J) 
[N=] 
IF(1.GV.KFC) IN=IMB—] 
AVOR=ABS(VOR(I,J)) 
YP=USTAR(II)*YNCIN,J)/VSC 
IF(YP.GT.50@.) THEN 
TI=YN(IN,J) 
ELSE 
T1=YN(IN,J)#(1.-EXP(-YP/26. )) 
ENDIF 
EDDY(1,J)=0.16*T1*T1*AVOR 
FY=AVOR «Tj 
IF(FY.GT.FMAX) THEN 
FMAX=FY 
YMAX=YN(IN, J) 
ENDIF 
ee ey a 
=@.5*(V(I,J+1)+V(1,J)) 
VSPHY=(UV*UV+VVeVV) /H(I, J) 
IF (VSPHY .GT .VDMX) VDOMX=VSPHY 
CONTINUE 
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C..OUTER EDDY VISCOSITY 


IF (FMAX.EQ.@.) GO TO 1@5 
FWAK E=F MAX * YMAX 
CFWAKE=@ . 25% YMAX # VDMX /FMAX 
IF(CFWAKE.LT.FWAKE) FWAKE=CFWAKE 
ICED=0 
DO 103 J=2,JL 
IF(II.LT. INOR1+1) be WORT) 
IF(II.GT.JP1) I=IMB-INOR(IMB—II, J) 
[N=] 
IF(I.GT.KFC) IN=IMB—I 
FKLEB=1./(1.4+5.50(@.3eYN(IN,J)/YMAX) #6. ) 
CEDDY=0 .01681.1*«FWAKE*sFKLEB 
IF(ICED.EQ.@ .AND. CEDDY.LT.EDDY(I,J)) THEN 
ICED=1 
IF(CEDDY/VSC.GT.@.) IOTUR(II)=1 
END I F 
IF(ICED.£Q.1) EDDY(1,J)=CEDDY 
103 CONT I NUE 
1805 CONTINUE 


C..MAKE SURE FLOW REMAINS TURBULENT AFTER TRANSITION 


1OT=@ 

IOTB=@ 

DO 110 I=KFC,1,-1 

1B=IMB—I 

reeecn 1M) IB=IM 

IF(IOTUR(I).£Q.1 .AND. IOT.£Q.@) I0T=I 
FAVOLENE 8) IOTUR(1I)=1 
IF(IOTUR(IB).EQ.1 .AND. IOTB.£Q.@) IOTB=IB 
IF(IOTB.NE.®) IOTUR(1B)=1 

11@ CONTINUE 


C..ASSIGN EDDY VISCOSITY IN THE WAKE ALONG WAKE GRID 


DO 120 JW=2, JR=-1 
IW=INOR(2, JW) 
1 BW=IMB—IW 
DO 119 J=JW,JUR 
I=IWK(JW, J) 
IF(J.EQ.INOR(2,J)) GO TO 116 
115 1B=]MB—] 
IF(I.£Q.1) IB=1 
EDDY( I , J)=EDDY(1W, JW) 
EDDY (1B, J)=EDDY( IBW, JW) 
116 I=I-1 
IF(IWK(JW-1,J).LT.1) GO TO 115 
119 CONTINUE 
120 CONTINUE 
RETURN 
END 


ibeps 


SUBROUTINE VORTY(IS,IL,INC, IC, ICBL,KODE) 


TeTe TET T TTT eS CST Tee ee TTT eee TET eee TTT Te TTT tee tee TT TTT Tt Te tt 
CALCULATE VORTICITY BY USING LINE- 
RELAXATION METHOD ON EACH RADIAL LINE 


CALLS OD TRID 


CALLED BY OD KNTCS 


RKSESHFESHSCSSHHSSCSC KK KK STK SSK SKK SEKKRKSKREKSKSESHKEKHSKRAKRERSEKKAREKSSEKESEHKSE DS 


MANDNDAAANNANA 


PARAMETER (IDIM=8@, JDIM=60) 
PARAMETER (KFC1=41) 


COMMON/COE/AF ,UI ,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 
COMMON/TUR1/USTAR(IDIM),EDDY(IDIM, JDIM) , IOTUR( IDIM), 

1 YN(KFC1,JDIM), IOT, IOTB 

COMMON/GRD/IM, IMZ,KFC,JR,N 
COMMON/VTEXT/KIN( IDIM) ,KEN(IDIM) 

COMMON /UNF /URR , URF , DFMX ,NCC 

COMMON/SGA/GA(IDIM, JDIM) ,GB(IDIM, JDIM) ,GC(IDIM, JDIM) 
COMMON/RHS/GD(IDIM, JDIM) ,DIP1(IDIM, JDIM) ,DIM1(IDIM, JDIM) 
COMMON/VV/VOR(IDIM, JDIM) , VOROLD( IDIM, JDIM) 
COMMON/SOLV/SUB(JDIM) , DIAG(JDIM) ,SUP(JDIM) , RHS(JDIM) 


IF(KODE.£Q.1) RETURN 
DO 10@ IC=1,NCC 
WMAX=0 . 
DMAX=@ . 
DO 11@ II=IS,IL,INC 
I=} ] 
IPPi=l+1 
IMizI—i 
Tae) 1M) I=] I—1M 


IF(1PP1.GT.IM) IPPi=IPPi-IM 
IF(IMi.GT.IM) IM1=sIMi-—IM 

JI=2 

IF(ICTUR.EQ.2 .AND. IOTUR(I).EQ.1) Jl=3 
JIMsJI—1 

JL=KIN(1) 


C..CONSTRUCT TRIDIAGONAL MATRIX AND SOLVE 


DO 120 JaJI,JL 
Ji=J-JIM 
SUP(J1)=GA(I, J) 
DIAG(J1)=GB(I,J) 
Sra Ne ee 
RHS(J1)=GD(1,J)+DIP1(1,J)*VOR(IPP1,J)+DIM1(I,J)*VOR(IM1,J) 
12@ CONTINUE 
RHS(1)=RHS(1)—-SUB(1)*VOR(I, JIM) 
CALL TRID(J1) 


C..UNDER-RELAX THE RESULT OF THE MATRIX SOLUTION 


DO 13@ J=JI,JL 

IF(ICBL .£Q. 1) THEN 
WW=RHS ( J—J IM) 

ELSE 
WW=URF *VOR(I , J )+URR®RHS(J—JIM) 
WABS=ABS (WW) 
DD=ABS (WW-VOR(I,J)) 
IF(DD.GE.DMAX) DMAX=DD 
IF(WABS.GT.WMAX) WMAX=WABS 


6 


ENDIF 
VOR(1,J)=Ww 
13@ CONTINUE 
118 CONTINUE 


C..IF BOUNDARY LAYER ZONE, RETURN TO KNTCS AFTER SOLVING EXPLICITLY. 
C ELSE ITERATE FOR CONVERGENCE AND RETURN WHEN CONVERGED OR WHEN 
C MAX NUMBER OF ITERATIONS IS EXCEEDED. 


IF(ICBL.E£Q.1) RETURN 
IF (DMAX/WMAX.LE.DFMX) RETURN 

10@ CONTINUE 
WRITE(6,10) I1,DMAX 
KODE=1 
RETURN 

10 FORMAT(2X,’°NO TRID CONVERGENCE AT I=’,14,’ DMAX=’,£10.3) 
END 


Loy 


QdONdAAARAIANON 


SUBROUTINE TRID(N) 


RESKSHREKHEKEKRHEKEERRHEKRREREKREREKREKEEKEEEREKEEKSEKEEEEEEREEEEEEREE EEE 


TRIDIAGONAL MATRIX SOLVER 
CALLS OD NONE 
CALLED BY OD VORTY 


SHEHEKCHKRHKHEKRSEKTKEKESHEEKSKEESEKHESSSSKEESSEESSERSEKSESEKREKEKESCHKCSEEHESESSEHKESEHEE 


PARAMETER (JDIM=60) 
COMMON/SOLV/SUB(JDIM) , DIAG(JDIM) , SUP( JDIM) ,RHS(JDIM) 


DO 100 J=2,N 
RATIO=-SUB(J)/DIAG(J—-1) 
DIAG( J )=DIAG(J)+RATIO#SUP(J—1) 
RHS (J )=RHS (J )+RATIO#RHS (J~—1) 
10@ CONTINUE 
RHS(N)=RHS(N)/DIAG(N) 
DO 110 J=N-1,1,—1 
11@ RHS(J)=(RHS(J)—-SUP (J) «RHS(J+1))/DIAG(J) 
RETURN 
END 


158 


OOAQDAAMAAAIAAAN 


SUBROUTINE FOCFT 


SRSESEHSSESeseeeSKeeCeeSKESESKseSsSesSseseesSeseseseeseseseseseseseseeseeseseseeeeeegues 


DETERMINE FOURIER COEFFICIENTS AT EACH RING 


CALLS @D NONE 


CALLED BY @D VORTY 


CPVAL 


SEEASSSS SEES SSSEESKCEESESESSSSSSESASSAa sess Seseseaese se see see sese ees 


180 


109 
118 


120 


PARAMETER (IDIM=8@, JDIM=60) 
PARAMETER (KFC1=41) 


COMMON/WFC/AA(JDIM, KFC1),BB(JDIM,KFC1) 
COMMON/SMT/SGMA(KFC1) 
COMMON/SCALE/H(IDIM, JDIM) 
COMMON/GRD/IM, IM2,KFC,JR,N 
COMMON/VTEXT/KIN(IDIM) ,KEN(IDIM) 
COMMON/TRIG/CS(KFC1, IDIM) ,SN(KFC1, IDIM) 
COMMON/FM/GAM( IDIM, JDIM) 


DO 100 J=1,N 
DO 10@ K=1,KFC 
BB(J,K)=0.0 
AA(J,K)=0.2 
CONTINUE 
DO 110 I=1,IM 
NI=KEN(1) 
DO 109 J=1,NI 
AA(J,1)@AA(J,1)+GAM(I,J) 
DO 109 K=2,KFC 
rN Sas UE corel ia baal ot 
BB(J,K)=8B(J,K)+GAM(I,J)*SN(K, I 
CONT INUE 
CONT INUE 
DO 120 J=2,N 
DO 120 K=2,KFC 
BAS Resa IK) a SeMALK) 
BB(J,K)=BB(J,K)*SGMA(K) 
CONTINUE 
RETURN 
END 


i) 


SUBROUTINE WSURF 


SESSTSESSETSEKSEKSSKSSESSK KK SKSEKRSEK SETHEKK SKEETER SKESSETESE 


COMPUTE SURFACE VORTICITY 
CALLS OD NONE 
CALLED BY @D KNTCS 


SSSSKKEKSHRHEKSEESHESTEHSCKSSSCCesseseSSesSSSKSSSseweeseseseSesee vee esses 


AQAAaAANAAAAANA 


PARAMETER (IDIM=8@, JDIM=60) 
PARAMETER (JP1=61) 
PARAMETER (KFC1=41,KFC2=42) 


COMMON/WFC/AA(JDIM,KFC1),88(JDIM,KFC1) 

COMMON /SMT /SGMA(KFC 1) 

COMMON/COE/AF ,UI,V1,OMG,VSC,NPL, ICTUR,ICST,ICPL 
COMMON/ABCDS/AS1 (KFC1) ,BS1(KFC1) ,CS1(KFC1) ,DS1(KFC1) 
COMMON/GRD/IM,IM2,KFC,JR,N 

COMMON/COR/RP(JP1,KFC2),RL(JP1) 

COMMON /VLB/VORLB ,NLB 

COMMON Q1(IDIM,JDIM),Q2(I1DIM,JDIM) ,Q3(IDIM, JOIM) ,Q4(IDIM,JOIM) 


C..DETERMINE FOURIER COEFFS. OF BOUNDARY VORTICITY, AA(1,M),BB(1,M) 
C FOR M.GE.4 


PIl=4. #ATAN(1. ) 
DO 100 K=4,KFC 
Cwi=0. 
DW1=0. 
DO 101 J=2,N 
ee OO dl a 
eet ieat | Kd) =W1 
OW1=0W1+8B(J ,K) =W1 
101 CONTINUE 
W2=1./(1./RP(1,K-3)—-1./RP(2,K=3) ) 
BB(1,K)=W2*((AS1(K )—DS1(K) )*FLOAT (K-3)-0W1) 
AA(1,K)=—W2*#((CS1(K)+8S1(K) ) *FLOAT (K-3)+CW1 ) 
10@ CONTINUE 


C..DETERMINE AA(1,K),88(1,K), FOR K.LE.3 


CW=0 . 

CW1=0. 

CW2=0. 

OW1=0. 

OW2=@. 

DO 130 J=2,N 
CW1=CW1+AA(J,2) #(RP(J+1,1)-RP(J,1)) 
OW1=0W1+8B(J ,2) «(RP(J+1,1)-RP(J,1)) 
CW2=CW2+AA(J,3)*#(RL(J+1)-RL(J)) 
DW2=D0W2+8B(J,3)*(RL(J+1)—-RL(J)) 

13@ CONTINUE 
DO 140 J=2,NLB 
140 CW=CW+AA(J,1)*(RP(J+1,2)—-RP(J,2)) 

AA(1,1)=—(4. #AF/P I *OMG+CW+2. *VORLB/P1)/(RP(2,2)—-RP(1,2)) 

AA(1,2)=(2.*VI-CS1(2)-BS1(2)—Cw1 )/(RP(2,1)-RP(1,1)) 

8B(1,2)=(AS1(2)-0S1(2)—-2. sUI-DW1 )/(RP(2,1)-RP(1,1)) 

AA(1,3)=—-(CS1(3)+BS1(3)+CW2)/RL(2) 

BB(1,3)=(AS1(3)—-DS1(3)-DW2)/RL(2) 

DO 15@ K=2,KFC 
AA(1,K)=AA(1,K) *#SGMA(K) 
BB(1,K)=BB(1,K) *SGMA(K) 

15@ CONTINUE 


160 





SUBROUTINE KMTCS 


LEE ER RE RE RE RERES SER EEE ERE RES ESE RSE R ESE ES EERE SE REESE SES ESE SE SE 


KINEMATICS OF THE PROBLEM 
CALLS OD NONE 
CALLED BY @0 ZONST 


SSHESHKSKSEHCHKSHLSSEKRKHKHEKCEKKSSSEEHAHEKRSEKSSEESHSSHEEHSKEKEEKEKEHESEKREKEKEESE 


QAAAAAAAIAN 


PARAMETER (IDIM=86, JDIM=60) 
PARAMETER (IP1=81, JP1=61) 
PARAMETER (KFC1=41,KFC2=42) 


DIMENSION VTOTPH(IDIM, JDIM) , VTOTCO(IDIM, JDIM) 


COMMON/COE/AF ,UI,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 
COMMON/VTEXT/KIN(IDIM) ,KEN(IDIM) 

COMMON/SMT /SGMA(KFC1 ) 
COMMON/ABCDS/AS1(KFC1) ,8S1(KFC1) ,CS1(KFC1) ,DS1(KFC1) 
COMMON/RHS/A(JP1,KFC1),B(JP1,KFC1),C(JP1,KFC1),D(JP1,KFC1) 
COMMON/WFC/AA(JDIM,KFC1) ,BB(JDIM,KFC1) 

COMMON /COR/RP(JP1,KFC2) ,RL(JP1) 
COMMON/VEH/UC(IP1,JP1),VC(IP1,JP1) 
COMMON/VV/VOR(IDIM, JDIM) , VOROLD( IDIM, JDIM) 
COMMON/VEL/U(IDIM, JP1),V(IDIM, JP1) ,HSTAR(IDIM, JDIM) 
COMMON/GRD/IM, IM2,KFC,JR,N 
COMMON/TRIG/CS(KFC1,1DIM) ,SN(KFC1, IDIM) 

COMMON Q1(IDIM,JDIM) ,Q2(IDIM, JDIM) ,Q3(IDIM, JDIM) ,Q4( IDIM, JDIM) 
COMPLEX HSTAR,VTOTCO,VTOTPH 


NP ah 1 
DO 106 K=i,KFC 
A rac K) 
B(1,K)=BS1(K) 
¢ hes! K) 
0(1,K)=DS1(K) 

1@@ CONTINUE 


DO 500 J=2,NP 
JM1=J—1 


C..BOUNDARY VELOCITY CONTRIBUTIONS 


Wi=.5/RP(J,2) 
C(J,2)=Wie2(C(1,2)—B(1,2) )+VI 
D(J,2)=Wi=#(D(1,2)+A(1,2))-UI 
A(J,1)=A(1,1)/RP(J,1) 
FAT OTR Re COU 
B(J,2)=Wi*(B(1,.2)-C(1,2))+VI 
DO 110 K=3,KFC 
Wi=.5/RP(J,K) 
C(J,K)=W1*(C(1,K)-B(1,K)) 
D(J,K )=W1*(D(1,K)+A(1,K)) 
A(J,K)=D(J,K) 
B(J,K)=-C(J,K) 
11@ CONTINUE 


C..INTEGRATE ALONG RADIAL AXIS FOR AA(J,K), BB(J.K); (K.GE.4) 
DO 120 K=4,KFC 
CW=0 


CW1=0. 
DW=O . 


G2 


DW1=0. 

DO 119 JJ=1,JM1 
W1=RP(JJ+1,K+1)—-RP( JJ, K+1) 
CW=AA(JJ,K) #W1+CW 
DW=8B( JJ ,K) «W1+DW 

119 CONTINUE 
IF(J.LT.NP) THEN 
DO 118 JJ=J,N 
W1=1./RP(JJ,K-3)—1. /RP(JJ+1,K-3) 
Sc HIE 
OW1=BB(JJ,K )*w1+DW1 
118 CONT INUE 

ENDIF 

W2=.5/(FLOAT(K+1)*RP(J,K)) 

CW2=W2 «CW 

DW2=W2+DW 

W3=.5eRP(J,K-2) /FLOAT(K—3) 

CW3=W3 «CW 

DW3=W3 «DW1 

Spas J.K)+CW2—CWS 

D(J,.K)=D(J,K)+DW2—Dw3 

A(J,K)=A Ha eee 

B(J,K)=B(J,K )-CW2-Cw3 

120 CONTINUE 


C..COMPUTE AA(J,K), BB(J,K):; (K.LE.3) 


CW=0. 

CW1=0. 

CW2=@. 

CW3=@ . 

CW4=0 . 

DWi=2. 

DW2=@ . 

DW3=8 . 

DW4=0 . 

DO 130 JJ=1,JM1 
CW=CW+AA (JJ, 1) #(RP(JJ+1,2)—RP( JJ ,2)) 
CW1=CW1+AA Tae eS are er ) 
OW1=D0W1+BB(JJ,2)*(RP TOs ce REN 
CW2=CW2+AA(JJ,3)e(RP(JJ+1,4)—-RP( JJ, 4) ) 
DW2=0W2+8B AUSSIE TS oRD eR 

13@ CONTINUE 

IF(J.LT.NP) THEN 

DO 146 JJ=J,N 
ea See re a rae ES Or 


OW3=0W3+8B( JJ ,2 ee JJ+1,1)-RP(JJ,1)) 
Picker 35. RL me ai 
DW4=0W4+8B (JJ.3)*(RL(JJ+1 )—-RL( JJ) ) 
140 CONTINUE 

ENDIF 

C(J,1)=(C(1,1)+.5eCw)/RP(J,1) 

CW3=.5*CW3 

DW3=.5*DW3 

W2=1./(6.*RP(J,2)) 

CW1=CW1 W2 

Dw1=DW1 »W2 


C(J,2)=C( J, 2)+CW1-CW3 
D(J,2)=0( J, 2)+DW1-DW3 
A(J,2)=A(J,2)+DW14DW3 
B(J,2)=8(J,2)—CW1-CW3 
w2=1./(8.*RP(J,3)) 
CW2=CW2 *W2 

DW2=DW2«wW2 


es 


W2=.5*eRP(J,1) 
CW4=CW4 # W2 

DW4=DW4 «W2 
C(J,3)=C(J,3)+CW2-CW4 
0(J,3)=0(J,3)+0W2-DW4 
A(J,3)=A(J,3)+DW2+0W4 
B(J,3)=8(J ,3)—CW2-CW4 

500 CONTINUE 


DO 15@ K=2,KFC 
DO 15@ J=2,N 
A(J,K)=A(J,K) #*SGMA(K) 
B(J,K)=B(J,K) #SGMA(K) 
C(J,K)=C(J,K)*SGMA(K) 
D(J,K)=D(J,K) #SGMA(K) 
15@ CONTINUE 


C..ABOVE COEFFICIENTS ARE FOR V(RHO) & V(PHI) IN NON-ROTATING 
C FRAME OF REFERENCE IN COMPUTATIONAL PLANE. 
C CALCULATE VELOCITIES WITH THESE VALUES. 


DO 16@ [=1,1M 
JL=KIN(1)+1 
IF(ICPL.£Q.0 .AND. NPL.GT.JL) JL=NPL 
IPP1=I+1 
IMi=I—1 
IF(1.£Q.1) IM1=IM 
IF(1.EQ.IM) IPP1=1 
Pa INURE ECE i) JL=KINCIPP1)+1 
IF(KINCIM1).GE.JL) JL=KIN(IM1)+1 
DO 159 J=2,JL 
ee OO 
V(I,J)*A(J,1)#.5 
DO 158 K=2,1M2 
ha Ese aed bar a ee ne eset eee 
V(I,J)=V(1,J)+A(d,K) CS(K,1)4B(J,K) &SN a 


158 CONTINUE 


Hen teat manatee Ue: = 
V(I,J)=V(1,d)+A(J,KFC) *CS(KFC,1)*.5 


159 CONT INUE 


C..TRANSFORM VELOCITY COMPONENTS TO ROTATING FRAME OF REFERENCE 
C FOR USE IN KNTCS 


DO 157 J=2,JL 
a aah gee pero 
W2=V(I ,J)*#SN(2,1)4+U(1,J)*CS(2,1) 
W1=W1+0MG#UC(1, J) 
W2=W2+0MG*VC(I1,J) 

IF( J.LE.60 nes 
VTOTCO(I,J) = CMPLX (W1,W2) 
VTOTPH(I,J) = VTOTCO(I,J)/HSTAR(I,J) 
Q2(1,J) = REAL(VTOTPH(1,J)) 

Q3(1,J) = AIMAG(VTOTPH(I,J)) 

ENDIF 

i Sched aa heap 

V(1,J)=W1*CS(2,1)+W2*SN(2,1) 


157 CONTINUE 
160 CONTINUE 


164 





DQANAIQAANANAANM| 


SUBROUTINE CPVAL 


SHESSESSSse—eeseeSseseseseseseaese_ esses seseseseseseeesvseesesese_seseseseseeeesease F 


CALCULATE CP VALUES BY CP INTEGRAL 


CALLS OD COED 


FOCFT 


CALLED BY @D ZONST 


SSSSHSHLSEESHCESSESKSESSSSSESCSSKSKRESSSKSSHSSSSEHESCHHESESLK TH ESSS 


10@ 


183 
105 


11 


Uae: 


PARAMETER (IDIM=80, JDIM=6@) 
PARAMETER (IP1=81 , JP1=61) 
PARAMETER Kroe 
PARAMETER (ICP=40) 


COMMON/GRD/IM, IM2,.KFC, JR,N 

COMMON /TRIG/CS(KFC1, IDIM) , SN(KFC1, IDIM) 

COMMON /VV/VOR(IDIM, JDIM) , VOROLD(IDIM, JDIM) 

COMMON/COE/AF,UI,VI,OMG,VSC,NPL, ICTUR, ICST, ICPL 

COMMON /SMT /SGMA (KFC 1 ) 

COMMON/VEL/U(IDIM,JP1),V(IDIM, JP1) 

COMMON /VTEXT/KIN(IDIM) .KEN(IDIM) 

COMMON/CPC/CCP(KFC1, JDIM) ,CP(1P1) 

COMMON/WFC/AA(JDIM,KFC1),BB(JDIM,KFC1) 

COMMON/RHS/A(JP1,KFC1),B(JP1,KFC1),C(JP1,KFC1) ,D(JP1,KFC1) 

COMMON/FM/GAM(IDIM, JDIM) 

COMMBN/CEM/ AW CTO 1B) BACT ER CCE) ONC I eG mace 
etic) ee ee ,CC2( ICP) ,DD2( ICP) 

COMMON Q1(IDIM, JDIM) ,Q2(IDIM,JDIM) ,Q3(IDIM, JDIM) ,Q4(IDIM, JDIM) 

DIMENSION re ctr aotoeal aaa r aa SLL A031) 
VCS(JDIM,KFC1),FC(KFC1) ,FS(KFC1) 


DO 105 I=1,IM 
JL=KEN(I) 
JLP=JL+1 
DO 100 J=1,JL 
GAM(I ,J)=VOR(1,J)/FLOAT(IM2) 
IF(JLP.LE.N) THEN 

DO 103 J=JLP,N 
GAM(I,J)=0. 

ENDIF 

CONT INUE 

DO 11@ J=1,KFC 
FC(1 )=@. 
FS(1)=@. 

CONT INUE 

CALL FOCFT 

DO 120 I=1,N 

I1= I+1 
DO 115 J=2,KFC 
Ji= J-1 
AW(J1)=AA(I,J) 
Bw(J1)=BB(I,J) 
CWw(J1)=.50e( C(I,J) + CC 
OW(J1)=.5e( D(I,J) + D(I1,J 
ath ee A(I,J) + ae 
FW(J1)=.59( B(I,J) +B 

CONT INUE 

AO=AA(I,1)#.5 

Com#25s( C(1,1) + €(11,1)) 

AOO=.25#( A(I,1) + A(I1,1) ) 

CALL COED(AO0,CO,A00,U0) 

UCC(I,1)=2. U0 


G6 


DO 116 J=2,KFC 
J1=J—1 
UCC(I, J)=AA2(J1) 
UCS(I, J)=8B2(J1) 
VEeECl ,J)ECE2( 31) 
vcS(I,J)=DD2(J1) 
116 CONTINUE 
120 CONTINUE 
DO 130 J=1,N 
DO 129 I=2,IM2 
Deets + rae mcrae en eCCP(I,J) 
FS(I)=FS(I) + (UCS(J,1)+VCC(J,1)) *CCP(I,J) 
129 CONTINUE 
FC(KFC)=FC(KFC) = UCC(J,KFC)*CCP(KFC, J) 
FC(1)=FC(1) + UCC(J,1)*CCP(1,J) 
130 CONTINUE 
DO 131 J=1,JDIM 
DO 131 I=1,IDIM 
Q4(1,J)=CCP(1,J)/.44.50(U(CI, J) ee24V(1,J0) 82) 
131 CONTINUE 
DO 140 I=2, IM2 
ne ecc(l) in eed 
FS(I)=FS(1) + VSC#AA(1,1 
140 CONT INUE 
DO 150 J=2,KFC 
. a ane 
FS(J)=FS(u)*SGMA(J) 
15@ CONTINUE 
DO 160 I=1,IM 
CP(I)=1.-FC(1) 
DO 159 L=2, IM2 
159 Ea eee 1) + 2.9( FC(L)*CS(L,1)=FS(L)*SN(L,1) ) 
CP(I)=CP(I) + FC(KFC)*CS(KFC,1) 
160 CONTINUE 
CP(IM+1)=CP(1) 
RETURN 
END 
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OOOINAANANQAAND 


SUBROUTINE COED(AO,CO,A00,U0) 


SSSSSEESSHEHSKSSKSHESHES SKEET SESSSSHSESEE SESE EKESEHEHSHEEEEEE EEE EEEEEEOE 


MULTIPLICATION OF TWO FOURIER SERIES 
CALLS OD NONE 
CALLED BY 0D CPVAL 


KREKSHSSSSESCSSS SSC KESKFESSSCHHSSCSKSKESKESKSESSTRHKAHTESSKEHSEKEHESKESEE 


PARAMETER (IDIM=80) 
PARAMETER (ICP=4@) 


COMMON/GRD/IM, IMZ,KFC,JR,N 
COMMON/CPW/AW(IDIM) ,BWCIDIM) ,CWC(ICP) ,DWCICP) , EW(ICP),FWC(ICP), 
1 AA2( ICP) ,BB2( ICP) ,CC2(ICP) ,0D2( ICP) 


NM1= IM2—1 
UO= COeA0 
DO 106 J=1,I1M2 
106 UO= UO + .58(AW( I) «CW(I)+6W(I)*OW(1)) 
DO 116 I=1,IM2 
II= [+I 
AA2(1I)= AW(II)*CW(I) + BWC(II)*DW(I 
BB2(I )= BW(II)*CW(I) —- ae ME 
CG Kz Se hea + BW(II)*FWCI 
002(1)= BW(II)*EW(1) ~ AW(II)*FW(I 
11@ CONTINUE 
DO 12@ J=1,NM1 
K4= IM2—J 
DO 119 K=i,J 
Kise IM2+i—K 
K2= Ki-K4 
K3= KitK4 
Wise AW(K2)+AW(K3 
W2= Bhike +BW(K3 
W3= AW(K2)—AW(K3 
W4= BW(K3)-BW(K2 
AA2(K4)= ar he + WisCW(K1) + MSeHiKaS 
BB2(K4)= BB2(K4) + W3eDW(K1) + W4eCW(Ki 
CC2(K4)= CC2(K4) + W1IeEW(K1) + ea 
DD2(K4)= OD2(K4) + W3*FW(K1) + W4sEW(K1 
11S CONT I NUE 
126 CONTINUE 
DO 13@ J=i ,NM1 
K4= J+1 
DO 129 K=1,J 
K2= K4—K 
K3= K4+K 
Wis Fs pee ee De 
W2= BW(K3)—-BW(K2 
W3= AW(K2)—AW(K3) 
W4= BW(K2)+BW(K3) 
AA2(K4)= AA2(K4) + W1IsCW(K) + Witenes 
BB2(K4)= BB2(K4) + W3sDW(K) + W4sCW(K 
CC2(K4)= CC2(K4) + WIieEW(K) + W2eFW(K) 
DD2(K4)= DD2(K4) + W3*FW(K) + W4eEW(K) 
129 CONTINUE 
13@ CONTINUE 
DO 148 I=i,IM2 
oe .5*AA2(1) + AOsCW(I) + eS ah 


BB2(1)= .5*BB2(1) + AO*DW(1) + CO*BW(1) 
CC2(1)= .5*CC2(1) + AO#EW(I) + ROREN CL) 
DD2(I)= .5*DD2(1) + AO#FW(1) + AOOsBW(I 
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14@ CONTINUE 
AA2(IM2)= 2. *AA2(IM2) 
CC2(IM2)= 2.*CC2(1M2) 
RETURN 
END 


bey, 


4e low la.O6o 4 
.0861 .0005 .0005 75 100 
Ace 
18 35 50 74 45 40 
.180 .10 50 
58 «650 O58 


ICST,RF,ALPS, ICTUR 

WMIN , DFMX , DRMX , KMAX , NCC 
URBI , URBP , URR 
IV1,181,182,1V2,NPL,NLB 
OT] ,DTINC ,NTMAX 
NTPL,NTOUT ,NTLO 
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MOOQAIAIAANIANAANAANAANAAAANDNN 


aQaaggqggNgn 


PROGRAM LOADS 


SSeeSeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeegeegeeguee 
LOADS COMPUTATIONS — DISSPLA VERSION 
LAST REVISION 4—7-87 


PRINCIPAL INVESTIGATOR : DR. J.C. WU 
AUTHORS : MIKE PATTERSON, ISHMAEL TUNCER 


GEORGIA INSTITUTE OF TECHNOLOGY 
(404) 894-3028 


TAPES : INPUT FROM GEOM 
TAPE4 : INPUT FROM ZONST 
TAPES : GENERAL INPUT 
TAPE6 : GENERAL OUTPUT 


TAPE1Q@ : OUTPUT TO PLOTS 
CALLS >: INTER 


SESHSHKSRLSESCKCRKSEKC ESE SCSKC SSCS KSC SHeSSSSSSSCSSesCe sees sSseseKC Se ese eseee eee s 


PARAMETER (IDIM=80, JDIM=60) 
PARAMETER (IP1=81, JP1=61) 


DIMENSION WB(IP1),CP(IP1),G(IP1) 
DIMENSION XB(IDIM),YB(IDIM) ,XT(IDIM) ,YT(IDIM) ,FUN(IP1) 
COMMON/L1/DTET, IM, IM1 


REWIND 3 
REWIND 4 
REWIND 5 
REWIND 6 
REWIND 10 


READ(5,*) NTOUT,NTS,NTL 
READ(3) AL,RE,DTET,XB,YB,XT,YT 


IM=IDOIM 
JR=JDIM 
IMi=IP1 
JR1=JP 1 
PI=4.*ATAN(1.) 


C..READ SURFACE VORTICITY & PRESSURE COEFFICIENT 


100 READ(4,END=99) NT,NTPL,T,RF,WB,CP,OMG,OMGD,ALP 


a eet NTS) GO TO 102 
IF(NT.GT.NTL) GO TO 99 
ALPD=ALP* 180. /PI 
WRITE(6,60) NT,T,ALPD 
pesecos (ALP) 
SINA=SIN(ALP ) 


C..COMPUTE TOTAL CP AND OUTPUT VALUES 


DO 110 I=1,I1M 


110 FUN(I)=OMGD* (XB(I)*YT(I)-YB(1I)*XT(I)) 


FUN(IM1)=FUN(1) 

CALL INTER(FUN,G,1,IM1) 

DO 12@ I=2,IM1 

WW=OMG * OMG* (XB(1) #XB(I)+YB(1)*YB(1)-XB(1)*XB(1)) 


12@ CP(1I)=CP(1)—G(1)-4W 


IF (MOD(NT,NTOUT).£Q.@) WRITE(6,61) (CP(I), I=1, IM) 


C..COMPUTE TANGENTIAL FORCE COEFFICIENT FROM VISCOUS EFFECT 


Lal 


DO 130 I=1,1M1 
13@ WB(I)=WB(1)—2.+*OMG 
DO 140 I=1,I1M 
140 FUN(1)=WwB(1)+«XT(1) 
FUN( IM1)=FUN(1) 
CALL INTER(FUN,G,1,IM1) 
CTV=2.@*G(IM1)/RE 


. COMPUTE NORMAL FORCE COEFFICIENT FROM VISCOUS EFFECT 


DO i856 I=1,IM 
15@ FUN(1)=WwB(1)*YT(I) 
FUN( IMi)=FUN(1) 
CALL INTER(FUN,G,1,1M1) 
CNV=2.0@G(IM1)/RE 


. COMPUTE MOMENT COEFFICIENT FROM VISCOUS EFFECT 


DO 160 I=i,IM 
16@ BUN (1a 8) Cy) 
FUN(IM1 )=FUN(1) 
CALL INTER(FUN,G,1,1M1) 
CMV=2. @eG(IM1)/(RE*AL) 
WRITE(6,66) CNV,CTV,CMV 


. COMPUTE TANGENTIAL FORCE COEFFICIENT FROM PRESSURE DISTRIBUTION 


DO 170 I=1,1IM 
17@ FUN(I)=CP(i)*YT(I) 

FUN( IM? )=CP(IM1)*¥T(1) 

CALL INTER(FUN,G,1,1M1) 

CTP=-G(IM1)/AL 


. .COMPUTE NORMAL FORCE COEFFICIENT FROM PRESSURE DISTRIBUTION 


DO 186 I=1,1M 
18@ FUN(I)=CP(1) eXT(1) 
FUN( IMi)=CP(IM1) #XT(1) 
CALL INTER(FUN,G,1,1M1) 
CNP=G(IM1)/AL 


. COMPUTE MOMENT COEFFICIENT FROM PRESSURE DISTRIBUTION 


DO 19@ I=1,IM 

19@ FUN(I)=CP(1)*(XB(1)*XxT(1)+YB(1)*YT(1)) 
FUN( IM1)=CP(IM1)¢(XB(1)*XT(1)+YB(1)*YT(1)) 
CALL INTER(FUN,G,1,1IM1) 
CMP=G(IM1)/(AL*AL) 


.. WRITE PRESSURE COEFFICIENTS AND COMPUTE TOTAL NORMAL, TANGENTIAL 
AND MOMENT COEFFICIENTS 


WRITE(6,64) CNP,CTP,CMP 
CN=CNV+CNP 
CT=CTV4CTP 
CM=CMV+CMP 


. COMPUTE LIFT AND DRAG COEFFICIENTS 
CLP=CNPsCOSA-CTPeSINA 
CDP=CNP«SINA+CTP*COSA 


CL=CNsCOSA-CT*SINA 
CD=CN*SINA+CT*COSA 
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WRITE(6,67) CLP,CDP,CMP 
WRITE(6,69) CL,CD,CM 


C..WRITE DATA FOR PLOTS PROGRAM AND RETURN TO BEGINNING OF PROGRAM 


WRITE(1@) NT,NTPL,T,ALPD,RF,RE,CLP,CDP,CMP,CP 
GO TO 108 


C..LAST OF DATA HAS BEEN READ 
99 WRITE(6,63) 











6@ FORMAT(1H ,40(1H-),/’ NT =',14,°, T =',F7.4 
+,’ AA =’, F8.4/1X,40(1H-)) 

61 FORMAT(//,’ CP VALUES, I=1,IM :’//(10E13.6)) 
63 FORMAT(//’? ——: END OF FILE :——'//) 
64 FORMAT(’ VALUES OF CNP,CTP,CMP ' 5X,3E16.8) 
66 FORMAT(//' VALUES OF CNV,CTV,CMV * 5X,3£16.8) 
67 FORMAT(’ PRESS. LOADS CLP,CDP,CMP —’ ,5X,3E16.8) 
69 FORMAT(’ TOTAL LOADS CL,CD,CM * 5X, 3E16.8) 

STOP 

END 


SUBROUTINE INTER(A,AI,IS, IL) 


SSSHSEHSHHASHASHHKESETSESEHRSHESHESESEAHEESHESEASEEHEEKEHTHTHERSESSHHEEEEKSEEE 


NUMERICAL QUADRATURE BY HIGH ORDER FORMULAS 
CALLS : NONE 
CALLED BY : LOADS 


SHESCKSeseSeeeessSeSeSCeesSKSEHSEKSHKSSSESSEHSSSSSKEKSSeeseseeeeeges 


MDAKOAAAAAQ 


DIMENSION A(IM1),AI(IM1) 
COMMON/L1/OTET, IM, IM1 


INC=1 
CI=DTET/24. 
IF(IL.LT.IS) THEN 
INC=~INC 
CIl=CI 
ENDIF 
AI(IS)=0. 
DO 10@ I=IS+INC,IL, INC 
IM2=I~2* INC 
IP1=I+INC 
IF(IM2.LT.1.OR.IM2.GT.IMi) IM2=ABS(IMi-1—IM2) 
IF(IP1.LT.1.0R.1P1.GT.IM1) IPi=ABS(IM1—1—IP1%) 
AI(1)=AI (I-INC)+CI*(—A(IM2)4+13. #(A(I-INC)4+A(1) )—A(IP1)) 
1@@ CONTINUE 
RETURN 
END 


ee 


PROGRAM PLOT1 


C SSEeeSeeeseeceeceeeeeeseegeeeeeeeeeseeeeeeeeeTCeeeecgegeeeeeengegeee 

C GEOMETRY PLOTTING PROGRAM 

C LAST REVISION 4-1-87 

é 

C PRINCIPAL INVESTIGATOR : DR. J.C. WU 

C AUTHORS : MIKE PATTERSON, ISHMAEL TUNCER 

C GEORGIA INSTITUTE OF TECHNOLOGY 

C (404) 894-3028 

C 

C  _TAPE1 : INPUT FROM GEOMETRY PROGRAM 

C TAPES : GENERAL INPUT 

C  TAPE11 : OUTPUT FOR PLOTTER 

C 

C CALLS : AXES, GRID, WAKE, ZONES 

C EQ KSKHSKHKESERKEKTKSSSSEOSKCSESCSSSESCSSCHKSESESHCHKESSEHSHOKSE SSE KSEE RVECKEKKEEHKECRE 
PARAMETER Be rhe) 
PARAMETER (IP1=81,JP1=61) 
PARAMETER (INOR1=20) 
COMMON /GRD/X(IDIM, JP1),Y(IDIM, JP1) 
COMMON/ZN/AL, JVL, JBL, JFL,1V1,1B1, IB2,1V2, ILINE(6) , SCALE 
COMMON/WK/INOR(INOR1, JDIM) , IWK(JDIM, JDIM) 
COMMON/PARAM/IM, JR 
DATA IRES/2/ 

c REWIND(1) 

C REWIND(S) 


READ(1) AL,X,Y,INOR, IWK 

READ(5,*) IV1,1B1,1B2,IV2 

READ sek IPLOTi , IPLOT2, IPLOT3, IPLOT4 
READ(5,«) JVL,JBL,JFL, SCALE, IDEV 


IM=IDIM 
JR=JDIM 


ILINE(2 )=IVi 
ILINE(3)=IB1 
ILINE s)o1v2 


Tce? 1 is 


ILINE(S)=#IV2 
ILINE(6)=IM 


C..SET DEVICE 
CALL DIP(15,’°PLOT1.DIP’ ,9) 


C IF(IDEV.EQ.0) THEN 

@ CALL EPIC(300. ,0,8.25,10.75,11) 
é ELSE IF(IDEV.EQ.1) THEN 
e CALL PCLCMP 

C CALL CALCMP(IBUF,512,11) 
C ELSE IF(IDEV.£Q.2) THEN 
¢ CALL PVRSTC 

fe CALL VRSTEC(IBUF,512,11) 
C ELSE IF(IDEV.EQ.3) THEN 
fe CALL P4115 

C CALL TK4115(IRES) 

C ENDIF 

C 


..CONSTRUCT EACH PLOT SPECIFIED 


174 


IF(IPLOT1.£Q.1) THEN 
CALL AXES(JVL,0.,’*VORTICITY GRID$’ ) 
CALL GRID(IM,IM,JVL) 
CALL ENDPL(@) 
END IF 
IF(IPLOT2.EQ.1) THEN 
CALL AXES(JBL,AL/4., AIRFOIL AND BL GRIDS’) 
CALL THKFRM(.@1) 
CALL FRAME 
CALL GRID(IM,IM,JBL) 
CALL ENDPL(@) 
END IF 
Baan ca 13 CALL ZONES 
IF(IPLOT4.EQ.1) CALL WAKE 


CALL DONEPL 

STOP 

END 

SUBROUTINE ZONES 


SSSSRSESSEKSKEHEKEKEEKSKEHHEEEHES KEES HKESESKEKEKEEKRSEKEKEEEKEEEEEEE 


AIRFOIL AND FLOW ZONES 
CALLS : AXES, GRID 
CALLED BY : PLOT1 


SSSSSSSeCeKeesesSeSseeeSeSeeK SSS SKS SKSHSHEHSEKEHEKREHESKEEEKKREKE SEES 


ONAQAMIAAAQAIAQO 


PARAMETER bane 
PARAMETER (1P1=81,JP1=61) 


COMMON/GRD/X(IDIM, JP1),Y(IDIM, JP1) 

COMMON/ZN/AL, JVL,JBL, JFL, 1V1,1B1, 1B2,1V2, ILINE(6) ,SCALE 
COMMON /PARAM/IM, JR 

DIMENSION XX(IP1),YY(IP1) 


CALL AXES(JFL,AL/4.,°FLOW ZONES$’) 
CALL THKFRM(.01) 

CALL FRAME 

CALL GRID(@,IM,5) 

CALL HEIGHT(.1) 


C..PLOT AND NUMBER THE DEMARCATION LINES 


DO 110 I=1,6 
DO 109 J=1,JFL 
a eX GILINE() 4} 
YY(J)=Y(ILINE(I),J) 
109 CONTINUE 
CALL pa ee Hite) 
CALL RLINT(ILINE(1),XX(JFL),Y(ILINE(I) ,JFL+1)) 
11@ CONTINUE 
CALL HEIGHT(.14) 
CALL ENDPL(@) 
RETURN 
END 
SUBROUTINE WAKE 


SSSSESSCKSCSSSSSSSSSSSTCSSSsSSSsSSSSSsSsesSseseseseseseseseseesee cece 


AIRFOIL AND WAKE GRID 
CALLS : AXES, GRID 


QAagagagang 


> 


C 
C 
C 


CALLED BY : PLOT1 


SRSHEKSKKESKREKHESEKESESSEKSEHEHSSHESESSEHSHEHEEHSSEHTHEEKREKSESESHEHEHEKEKESEEEEE 


PARAMETER (IDIM=8@, JDIM=6@) 
PARAMETER (IP1=81,JP1=61) 
PARAMETER (INOR1=20) 


COMMON/GRD/X(IDIM, JP1),Y(IDIM, JP1) 

COMMON/ZN/AL, JVL, JBL, JFL, 1V1,1B1, 1B2, 1V2, ILINE(6) , SCALE 
COMMON/WK /INOR(INOR1,JDIM) , IWK(JDIM, JDIM) 
COMMON/PARAM/IM, JR 

DIMENSION XX(IP1),YY(IP1) 


C..PLOT AXES AND AIRFOIL 


CALL AXES(JVL,AL, ’WAKE GRID$’ ) 
CALL GRID(@,IM,5) 


C..PLOT VORTICITY GRIO UNDERNEATH WAKE GRID 


CALL DOT 

CALL THKCRV(.@1) 

CALL GRID(IM/4+1, IM/4+1, JR) 
CALL RESET(’DOT’) 

CALL RESET( ’ THKCRV’ ) 


C..PLOT STREAMWISE LINES 


DO 11@ JN=1,JUR~1 

JJ=O 

DO 109 J=JN+1,JR 
JJ=JU+i 
a TSeeL ONC 
YY (JJ )=YCIWK(JN,J) J) 


CALL CURVE(XX, YY, Ju,@) 


1@8 CONTINUE 


11@ CONTINUE 


C..PLOT NORMAL LINES 


NOOIOAIAIAANAYN 


DO 126 J=2,IM/4 

DO 119 J=1,JR 
YE errpGetn: 34.35 
YY(J)=YCINOR(I,J),J) 


CALL CURVE(XX, YY, JR,@) 


119 CONTINUE 
120 CONTINUE 


CALL ENDPL(@) 

RETURN 

END 

SUBROUTINE AXES ( JMAX, OFFSET, LABEL) 


SBRSESSKESESHSHSKEKEHSEHHESEHKEKESKREKESEKEESHETEKESEKSSECHEHKEKERESECHKEER 


DETERMINE THE LAYOUT OF THE PAGE 
CALLS : NONE 


CALLED BY : PLOT1, ZONES, WAKE 


SHHSSKSHKHESSEHKSEHEKRSEHSHSEKSHK SESE HEKESSKCHTSEHKEKLSHHHKKEHEESSSKEKSKSTKBEKSETS 


PARAMETER (IDIM=8@, JDIM=6Q) 
PARAMETER (IP1=81,JP1=61) 


COMMON/GRD/X(IDIM, JP1),Y(IDIM, JP1) 
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Cc. 


Cc. 


C 


Cc. 


QAMANDANAND 


COMMON/ZN/AL, JVL, JBL, JFL, 1V1, 161, 182, 1V2, ILINE(6), SCALE 


COMMON /PARAM/IM, JR 
CHARACTER*2@ LABEL 


.PAGE LAYOUT 


PAGEX=8 . 25 

PAGEY=10.75 

XAXIS=7.25 

YAXIS=7.25 

CALL BLOWUP (SCALE) 

CALL PAGE(PAGEX,PAGEY) 

CALL AREA2D(XAXIS, YAXIS) 

CALL HEADIN(LABEL, 100,1.25,1) 


.COMPUTE GRID AXIS LENGTHS AND SCALES BASED ON 
MAXIMUM VALUES TO BE PLOTTED 


GRIDMAX=0. 
DO 116 I=1,IM 
DO 109 J=1,JMAX 
eat .GT. GRIDMAX) GRIDMAX=ABS(X(I, J 
IF(ABS(Y(I,J)) .GT. GRIDMAX) Sa ReaBaiee eas 
189 CONTINUE 
11@ CONTINUE 


. LOCATE THE ORIGIN OF THE PLOT AND DRAW AXES 


XOR IG=—GR I DMAX+0F FSET 

YOR 1G=--GR I DMAX 

XMAX=GR I DMAX+OF FSET 

YMA X=GR I DMAX 

XSTP=1. 

YSTP=1. 

CALL GRAF (XORIG,XSTP,XMAX, YORIG, YSTP, YMAX) 
RETURN 

END 

SUBROUTINE GRID(IMAX1, IMAX2, JMAX) 


SBSESSeKCSKRESKKESKKESSKSSSKESKSEKSEKCKEKESKHSSEKR SESS SKSEKSESSEeSE KEE 


GRID IN PHYSICAL PLANE 
CALLS : NONE 
CALLED BY : PLOT1, ZONES, WAKE 


SKSSTKKSERHSSEKSKEKSKSKKSKHSSEHSKEEKSEKKSEKHHEHRSHEKEKEKKKESESEE KEKE LE 


PARAMETER (IDIM=8@, JDIM=60) 
PARAMETER (IP1=81,JP1=61) 


COMMON/GRD/X(IDIM, JP1), Y(IDIM, JP1) 
COMMON/PARAM/IM, JR 
DIMENSION XX(IP1),YY(IP1) 


.PLOT THE RADIAL LINES 


DO 110 I=1, IMAX1 
DO 109 J=1, JMAX 
XX(J)=X(I,J) 
YY(J)=Y(I,J) 
109 CONTINUE 
CALL CURVE(XX,YY, JMAX,@) 
11@ CONTINUE 


IL yi 


C..PLOT AZIMUTHAL LINES 


I = IMAX2 
DO 12@ J=1,JMAX 
DO 119 I=1, IMAX2 
XX(I)=X(1,J) 
YY(1)=Y(1,J) 
119 CONTINUE 
IF (IMAX2.EQ.IM) THEN 
le i eae a 
YY (IMAX2+1 )=YY(1) 
I 1=]MAX2-+1 
END IF 
CALL CURVE(XX,YY, 11,0) 
12@ CONTINUE 
RETURN 
END 


ses 


ONNAANAANANAANAANANANANAANAN 


QaQaq0 


Oaoaana 


PROGRAM PLOT2 


SSEASKKSEKASASKaAaSsASseKSeSHAAHAAHSAHSHSSAHAsAsAsSsasceSs sas HAeacseeseeassaasaasgage 
THIS PROGRAM ACCEPTS 2D-ARRAYS TO MAKE CONTOURS 
LAST REVISION 5-25-86 
PRINCIPAL INVESTIGATOR : DR. J.C. WU 
AUTHOR : MIKE PATTERSON 
GEORGIA INSTITUTE OF TECHNOLOGY 
(404) 894-3028 
TAPES : GENERAL INPUT 
TAPE6 =: GENERAL OUTPUT 
TAPE9 =: INPUT FROM ZONST FOR STREAMLINES AND CONTOURS 


TAPE12 : OUTPUT FOR PLOTTER 
TAPE14 : INPUT FROM GEOM 


CALLS : CONTG 


SKSASSSASCHASASHAKASSHSHKHKESHASEKHKHKHSSAAHSLSEAAASSSEKSKAHKEKEKRKKEKSKE RHE 


PARAMETER (IDIM=8@, JDIM=60) 
PARAMETER (IP1=81, JP1=61) 


REAL XX(IP1,JP1),YY(IP1,JP1),S1(IP1,JP1),S2(IP1,UP1) 
REAL XB(IDIM),YB(IDIM) ,CONT(1@0) 

REAL XB1(IP1),YB1(IP1) 

INTEGER IBUF(512) 
COMMON/SUM/CSQ, GAMA, SIGMA, COSA, SINA 

DATA IRES/2/ 

DATA PI/3.14159/ 


REWIND 5 
REWIND 6 
REWIND 9 
REWIND 14 


CALL DIP( 15, ’PLOT2.DIP’, 9 ) 
.. GENERAL INPUT 


READ(5,*) VALMIN, VALMAX ,NCON 

READ(5,*) VORMIN, VORMAX ,NVOR 

READ(5,*) AHI, ICHR,NDIG,NSKIP, LINTP, HEAVY 
READ(5,*) X1,X2,DX,XLEN 

READ(5,*) Y1,Y2,DY,YLEN 

READ(5,*) IDEV, IOPT, IPRINT,IPAR 

READ(5,*) JMAX,SCALE,NTS,NTL 
IF(NCON.GT.@) READ(5,*) (CONT(I),I=1,NCON) 


IM=IDIM 
IM1=IP1 


.. INPUT FROM GEOM 


READ(14) CSQ,GAMA,SIGMA,AL 
READ(14) XB, YB,XX,YY 


foe! DEVICE 


IF(IDEV.EQ.@) THEN 
CALL EPIC(30@. ,0,8.25,10.75, 12) 
ELSE IF(IDEV.EQ.1) THEN 
CALL PCLCMP 
CALL CALCMP(IBUF ,512, 12) 
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ELSE IF(IDEV.£Q.2) THEN 
CALL PVRSTC 
CALL VRSTEC(IBUF, 512,12) 
ELSE IF(IDEV.£Q.3) THEN 
CALL P4115 
CALL TK4115(IRES) 
ENDIF 

CALL SETDEV(6,6) 


aQgagaaqnnnaan 


C..INPUT FROM ZONST —= OUTER LOOP — NEW PAGE OF PLOTS 


1@0@ CONTINUE 

READ(9,END=3000) NT,NPL,T,ALPD,RF,RE,S2,S1 

IF(NT.LT.NTS) GO TO 1000 

IF(NT.GT.NTL) GO TO 3000 

WRITE(6,1@) NT,NPL,T,ALPD,RF,RE 

IF(JMAX.GT.NPL) THEN 
WRITE(6,%) ‘JMAX > NPL ——EXPECT GARBAGE——’ 
STOP 

ENDIF 

ALP=ALPDsP1/18@. 

ab eT. 

COSA=COS (ALP 

1SUB=@ 

INAME=IOPT 

CALL NOBRDR 

CALL NOCHEK 

CALL GRACE(@. ) 


C..INNER LOOP — SUBPLOTS 


200@ CONTINUE 
1SUB=1SUB+ 1 
fa sey CALL PHYSOR(1.5,5.50) 
IF(1SUB.EQ.2) CALL OREL(@. ,-4.25) 
CALL COMPLX 
CALL BASALF(’STAND’ ) 
CALL MIXALF(’L/CSTD’ ) 
CALL MX3ALF(’L/CGR® , 1H/) 


CALL BLOWUP(SCALE) 
. CALL AREA20(XLEN, YLEN) 
CALL YAXANG(@.@) 
CALL XTICKS(@) 
CALL YTICKS(@) 


CALL XNONUM 
CALL YNONUM 
CALL GRAF (X1,DX,X2,Y1,DY, Y2) 


CALL THKFRM(.@2) 
CALL FRAME 


C..PLOT PARAMETERS 


IF(IPAR.EQ.1) THEN 


WB = 1.55 
HB = 1.@5 
XBLK = .@5 
YBLK = .@5 


CALL BLTREC(XBLK, YBLK,WB,HB,@.0,0.215) 
CALL BECEM(ID) 
CALL BLOFF(ID 
XM = XBLK + 0.10 


ire 0 


YM = YBLK 


+ 0.05 


CALL HEIGHT(@. 15) 
CALL MESSAG(’R(E =)$',100,XM, YM) 


XR=XM+ .5 


CALL REALNO(RE,-1,XR, YM) 


XM = XBLK 
YM = YBLK 


orale a 
a0 250 


CALL HEIGHT(Q. 15) 

CALL MESSAG(’(RF =)$°,100,XM, YM) 
XR = XM + 0.6 

CALL REALNO(RF,3,XR, YM) 


XM = XBLK 
= YBLK 


0.1 
te oo) 


CALL HEIGHT(@. 15) 

CALL MESSAG(’(T = =)$’°,100,XM, YM) 
XR = XM + 0.6 

CALL REALNO(T,3,XR, YM) 


XM = XBLK 
YM = YBLK 


+ 0.1 
+ 0.80 


CALL HEIGHT(@. 15) 

CALL MESSAG(*/A)  =$°,10@,XM, YM) 
XR = XM + 0.6 

CALL REALNO(ALPD,3,XR, YM) 


CALL BLON 


(1D) 


C..ORAW BLADE AT ANGLE OF ATTACK 


140 


DO 140 I=1, 
XB1(1)=XB 
YB1(1)=YB 

CONT INUE 


IM 
eae, WEE Bets 
I) *COSA-XB(I)*SINA 


XB1( 1M1)=XB1(1) 


YB1(IM1)=YB 
CALL THKCRV 


1(1) 
(8.010) 


CALL CURVE(XB1,YB1,IM1,Q) 


CALL RESET( 


* THKCRV’ ) 


C..DISPLAY DATA, CALL CONTOURING ROUTINE 


i 
2 


1 


1 


IF(IPRINT.EQ.1) 


WRITE(6,20)((1.u,XX(I,d).YY(I,d 


IF(ISUB. EQ. 


CALL CONTG(S1,XX,YY,1,1M1,1, JMAX, VALMIN, VALMAX, AHI ,CONT,NCON 


EESE 


J=i,JMAX,2),1=1,1M,4 
1) THEN 
ICHR ,NDIG,NSKIP, LINTP , HEAVY) 


genes aCe) 


CALL CONTG(S2,XX,YY,1,1M1,1,JMAX, VORMIN, VORMAX , AHI ,CONT,NVOR, 


ENDIF 


ICHR ,NDIG,NSKIP,LINTP, HEAVY ) 


C..WRITE APPROPRIATE JITLE 


XM=XLEN—4 .@ 
YM=YLEN-.25 
IF(INAME.£Q.-2) THEN 


CALL MESSAG(’ S(TREAMLINES)$’, 100, XM, YM) 


ENDIF 
IF(INAME.EQ.-1) THEN 


CALL MESSAG(’V(ORTICITY) C(ONTOURS)$’ , 100,XM, YM) 


ire 


ENDIF 
IF(INAME.EQ.@) THEN 

CALL MESSAG(’D(ENSITY) C(ONTOURS)$’ , 100, XM, YM) 
ENDIF 
IF(INAME.EQ.1t) THEN 

CALL MESSAG(’M(ACH) N(UMBER) C(ONTOURS)$’ , 180, XM, YM) 
ENDIF 
IF( INAME.EQ.2) THEN 

CALL MESSAG('V(ELOCITY) M(AGNITUDE) C(ONTOURS)$° , 102, XM, YM) 
ENDIF 
IF(INAME.EQ.3) THEN 

CALL MESSAG(’S(KIN) F(RICTION) C(ONTOURS)$’ , 100, XM, YM) 
ENDIF 


C..FINISH SUBPLOTS OR START NEW PAGE 


CALL ENDGR(ISUB) 

J NAME= I NAME-+ 1 

IFC ISUB.EQ.1) GO TO 2006 
CALL ENDPL(@) 

GO TO 1060 


C..PLOTTING FINISHED 


NAXAQAANQNANDAN|AOA 


30@6 CALL DONEPL 


STOP 
1@ eet ’, NT,NPL,T,ALPD,RF,RE’ ,214,4F13.3) 
2@ FORMAT(2(5X,215,4F1@.5) ) 

END 


SUBROUTINE CONTG(A,XA,YA,IL,1U,JL,JU,VALMIN, VALMAX, 
AHI ,CONT,NCON, ICHR, NDIG,NSKIP,LINTP, HEAVY) 


SERHSOKSOKSLHSSSOHSSHH SSSR HTHESSHHTHE SESS HSE SEE AHES SESE SEEEEE EE 
THIS SUBROUTINE PLOTS CONTOURS OF MATRIX "A" FOR GENERAL (X,Y) 
POSITIONS OF THE ELEMENTS 


CALLS : XFUN 


CALLED BY : PLOT2 


SSSSCESKCESESEKSSSSSSKSESCSSSHSCKTSKEESHSESHSSCSSSSE SS esesSsseseseseseseeseseseseasesee ve 


PARAMETER (IP1=81, JP1=61) 

DIMENSION AIDE Gaia C0 loc 28. 1s) 
DIMENSION B Pi ae 

COMMON /COORD/XX (3), YY(3) ,XX2(3) ,YY2(3) 

REAL HT(4) 

LOGICAL HVY 

DATA HT/.08,.16,.12,.14/ 

DATA PI/3.14159/ 


IF(ICHR.GT.@) CALL HEIGHT (HT(ICHR) ) 

12=IU-1 

J2=JU-1 

NSKP=NSKIP 

IF(LINTP.EQ.3) CALL CHNDOT 

IF(LINTP.EQ.4) CALL CHNDSH 

IF(LINTP.EQ.2) CALL DASH 

IF(LINTP.£Q.1) CALL DOT 

IF(NSKIP.LE.@) NSKP=10 

HVY=. FALSE. 

IF(VALMIN.LT.VALMAX) THEN 
VALMX=VALMAX 
VALMN=VALMIN 

ELSE 
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VALMN=1£30 
VALMX=-1£3@ 
DO 170 I=IL,IU 
DO 170 J=JL, JU 
IF (A(I,J).LT.AHI) THEN 
VALMX=AMAX1(A(I,J) , VALMX) 
VALMN=AMIN1(A(I, J), VALMN) 
ENDIF 
170 CONT INUE 
ENDIF 
N=NCON 
IF(NCON.LT.@) THEN 
= —-NCON 
CONMIN = VALMIN 
CONINC = (VALMAX — VALMIN)/FLOAT(N-1) 
ENDIF 


C..DO LOOP FOR CONTOUR LINES 


DO 5@@ NCN=1,N 
IF(NCON.LT.@) THEN 
CONTUR=CONMIN+(NCN—-1 ) #*CONINC 
ELSE 
CONTUR=CONT (NCN) 
ENDIF 
IF (HEAVY .NE.@. ) 
1 HVY=ABS(AINT (ABS (CONTUR/HEAVY )+0 .5)—-ABS (CONTUR/HEAVY) ) .LT.@.2@1 
IF(HVY) CALL THKCRV(.01) 
NSK=NSKP 


C..SEARCH MATRIX "A" FOR CONTOURS 


DO 490 I=IL,12 
DO 490 J=aJL,J2 
BB=AMAX1(A(1I,J),AC1,J+1) ,ACI+1,J0) ,ACI+1, J#1)) 
IF (BB-AHI) 290,492,490 
290 IF (CONTUR-BB) 300, 490,490 
30@ IF (CONTUR-AMIN1(A(I,J),A(I,J+1),A(1+1,J),AC(I+1,0+1))) 490,310,310 
310 B(1)=.25#(AC(I,J)+AC1I+1,0)+AC1, J+1)+A( 141, +1) ) 
B(4)=B(1) 
X(1)=.252(XACI , J )+XAC 141, 5)+XAC 1, S41 )+XA( 141, J+1)) 
X(4)=X(1) 
Y¥(1)=.258(YACI, J)+YACI4+1,J5)+YACI , J+1 )4+VAC 141, J+1)) 
Y(4)=Y(1) 


C..SEARCH MATRIX SUB-CELL TRIANGLES FOR CONTOURS 


DO 48@ K=1,4 
NP=@ 

IF (K.LE.1) THEN 
B(2)=A(I+1,J) 
X(2)=XA(I+1,J) 
Y(2)=YA(I+1,J) 
E 


B(2)=B(3) 
tee 
¥(2)=¥(3 
ENDIF 
IM=I+K/3 
JM=J+1—IABS (5—2%K) /3 
8(3)=A( IM, JM) 
SO Ua 
Y(3)=YACIM, JM 


Jlrs 


C..DETERMINE CONTOUR INTERSECTIONS 


DO 430 M=1,3 
IF (CONTUR-AMIN1 (B(M).B(M+1))) 430,380, 380 
380 IF (CONTUR-AMAX1(B(M),B(M+1))) 390,390, 430 
396 NP=NP+1 
BB=B(M+1)—B(M) 
IF (ABS(BB).LE.1.€-15) THEN 
D=@.5 
ELSE 
D=(CONTUR-8(M) ) /BB 
ENDIF 
SE eae ea 
YY (NP) =Y (M)+D*(Y (Mei )-Y(M) ) 
43@ CONTINUE 
IF(NP.GT.1) THEN 
IF(ICHR.LT.@) THEN 
ba (oo ae ae CALL RESET(’DASH' ) 
IF(CONTUR.LE.@.) CALL DASH 
ENDIF 


C..TRANSFORM CONTOUR TO PRYSICAL PLANE AND DRAW 


CALL XFUN(NP) 
CALL CURVE(XX2, YY2,NP,@) 


C..LABEL CONTOURS 


IF(ICHR.GT.@) THEN 
NSK=NSK+1 
IF(NSK.GE.NSKP) THEN 
NSK==1 
XS=. 5 (XX2(1)+XX2(2)) 
DXS=XX2(2)~XX2(1) 
YS=.5*(YY2(1)+YY2(2)) 
DYS#YY2(2)—YY2(i) 
ANG=18@. *ATAN2(DXS,DYS) /PI 
CALL ANGLE(ANG) 
CALL REALNO(CONTUR,NDIG,XS,YS) 
ENDIF 
ENDIF 
ENDIF 
48@ CONTINUE 
49@ CONTINUE 
IF(HEAVY .NE.@) CALL RESET(’ THKCRV’ ) 
5@@ CONTINUE 
RETURN 
END 
SUBROUTINE XFUN(NP) 


CHHKSHSHAKEH SHEESH KRESS SESE SHEESH ASR SESE STEER ERESES 
COMPUTES CARTISIAN COORDINATES (XX,YY) IN PHYSICAL PLANE 
OF A SPECIFIC POINT (X,Y) IN COMPUTATIONAL PLANE 


CALLS >: NONE 
CALLED BY : PLOT2, CONTG 


SSSSeeseseseseses Sees se Cee ss sess SC seseses SC Sess sess eseeseeseseseseese cesses 


MNAAMQAAIIAAAN 


COMMON /SUM/CSQ , GAMA, SIGMA,COSA,SINA 
COMMON/COORD/X(3),Y(3) ,XX(3),Y¥(3) 
COMPLEX Z1,Z 


DO 10@ I=1,NP 
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Z1=GAMA+CMPLX(X(1),Y(1)) 
Z=Z1+CSQ/Z1+SIGMA 
X1=REAL(Z) 
Y1=AIMAG(Z) 
XX (1 }=X1*COSA+Y1*SINA 
YY (1 )=Y1*#COSA-X1*SINA 
100 CONTINUE 
RETURN 
END 


ESS 


MQQAQANDADAIAIAI2NRI|A|AQIAAAIANA 


NO ONQQNdNgNdNNdNd1NND OY 


PROGRAM PLOTS 


SHAKER AAAS TAA TERRA RRR TERRE RR RRR ETE 
LOADS PLOTTING PROGRAM 
LAST REVISION 6-25-87 
PRINCIPAL INVESTIGATOR : DR. J.C. WU 
AUTHOR : MIKE PATTERSON 
GEORGIA INSTITUTE OF TECHNOLOGY 
(464) 894-3028 
TAPES : GENERAL INPUT 
TAPE6 > GENERAL OUTPUT 


TAPE1@ : INPUT FROM LOADS PROGRAM 
TAPE13 : OUTPUT TO PLOTTER 
TAPE14 : INPUT FROM GEOM 


CALLS : NONE 


SESE KSKESHEHKSESTEPSHKRSESESSESLHSKKEKSSKSRESHEFHSSCCRESESECESCBERE 


PARAMETER (JDIM=8@, JDIM=60) 
PARAMETER (IP1=81,JPi1=61) 
PARAMETER ([U=40, IL=42, IMAX=200@ ) 


REAL XB(IDIM) ,XU(IU) .XL(IL) 

REAL CP(IP1),CPU(IU) ,CPL(IL) 

REAL T(IMAX),CL(IMAX) ,CD(IMAX) ,CM(IMAX) , ANG( IMAX) 
REAL XEXP(IU), CPEXP(IU,JDIM), CLEXP(IU), CDEXP(IU), 
> CMEXP(IU), ALPHA(IU) 

INTEGER IBUF(512) ,NT(IMAX) 

DATA IRES/2/ 


REWIND 5 
REWIND 6 
REWIND 16 
REWIND 14 


. READ INPUT PARAMETERS 


READ(5,*) SCALE,SCALE2 

ROS X1,X2,DX,XLEN 

READ(S,*) Yi,Y2,DY,YLEN 
READ(5,*) IOPT,IDEV, IPAR,NTS,NTL 


.. INPUT FROM GEOM 


READ( 14) DUMM1,DUMM2,DUMM3,AL 
READ(14) xB 


.. INPUT FROM EXPERIMENTAL RESULTS 


READ(24,*) NPARTS 
READ(24,900) ( XEXP(1),1=1,16 ) 
DO J=1,NPARTS 
READ(24,90@) ALPHA(J). CLEXP(J), CDEXP(J), CMEXP(J) 
READ(24,90@) ( CPEXP(I,J),1=1,16 ) 
DO I=1,16 
CPEXP(I,J) = —CPEXP(I,J) 
ENDDO 
ENDDO 


IEXP = @ 
900 FORMAT( 5( 1X,£14.7 ) ) 
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C..SET DEVICE 
CALL DIP( 15, ‘PLOT3.DIP’, 9 ) 


IF(IDEV.EQ.@) THEN 
CALL EPIC(300. ,0,8.25,10.75,13) 

ELSE IF(IDEV.EQ.1) THEN 
CALL PCLCMP 
CALL CALCMP( IBUF ,512, 13) 
ELSE IF(IDEV.£Q.2) THEN 
CALL PVRSTC 
CALL VRSTEC(IBUF ,512, 13) 
ELSE IF(IDEV.£Q.3) THEN 
CALL P4115 
CALL TK4115(IRES) 
ENDIF 

CALL SETDEV(6,6) 

IM=IDIM 


DQANQNANAANAIANNO 


C saseassssCP PLOT LOOPssasccccsccs 
IF(IOPT.£Q.@) THEN 
C..INPUT FROM LOADS 
10@@ CONTINUE 
READ(10,END=200@) NT1,NTPL,T1,ALPD,RF,RE,CL1,CD1,CM1,CP 


IF(NT1.LT.NTS) GO TO 1000 
IF(NT1.GT.NTL) GO TO 2000 


c IF (MOD(NT1,NTPL).NE.@) GO TO 1000 
WRITE(6,1@) NT1,71,ALPD 
1@ FORMAT(’ °,’NT,T,ALPD :’,15,2F8.3) 
Cc IEXP = IEXP + 1 


C..SET UP ALPHA-NUMERICS AND AXES NAMES 


CALL NOBRDR 
CALL COMPLX 

CALL BASALF(’STAND’ ) 
CALL MIXALF('L/CSTD’ ) 
CALL MX3ALF(’L/CGR’ , 1H\) 
CALL BLOWUP(SCALE) 


C..LABEL AXES 


CALL qe, 178 128) 
CALL YNAME(’-C(P)$’ , 100) 


C..DIVIDE DATA INTO UPPER AND LOWER ARRAYS 


DO 10@ I=1,IU 
XU (I )=(XB(I)+AL/4.)/AL 
CPU(I)=CP(I) 
10@@ CONTINUE 


C IF( ITEST.£Q.@ )THEN 

c WRITE(24,*) ( XU(I),I=IU,1,-1 ) 
C ITEST=1 

C ENDIF 


WRITE(24,*) ALPD 
WRITE(24,*) ( CPU(I),I=IU,1,-1 ) 


IM2=IM/2—1 
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DO 110 I=1, IL-1 
XL(1)=(XB(1+IM2)+AL/4.)/AL 
CPL(1)=-CP(I+IM2) 

110 CONTINUE 

XL(IL)=xU(1) 

CPL(IL)=CPU(1) 


C..DEFINE PLOTTING AREA AND HEADING 


CALL PHYSOR(1.75,3.0) 

CALL AREA2D(XLEN, YLEN) 

CALL SCLPIC(SCALE2) 

CALL YAXANG(0.@) 

CALL XTICKS (5) 

CALL YTICKS (5) 

CALL GRAF (X1,DX,X2,Y1,DY,Y2) 
CALL THKFRM(.@2) 

CALL FRAME 

CALL RESET(’ THKFRM’ ) 


C..PARAMETERS 


WB = 2.65 

HB = .25 

XBLK = XLEN-3.5 

YBLK = YLEN+O.4 

CALL BLTREC(XBLK,YBLK,WB,HB,@.0,0.015) 
CALL BLKEY(ID1) 

CALL BLOFF(ID1) 

XM = XBLK + 0.10 

YM = YBLK + 0.05 

fe CALL peLenre ts) 

C CALL MESSAG(’P(RESSURE) D(ISTRIBUTION)$‘ ,1@@, XM, YM) 


aan 


WB = 1.55 
HB = 1.05 
XBLK = XLEN-1.6 

YBLK = YLEN-1.1 

CALL BLTREC(XBLK, YBLK,WB,HB,@.0,0.015) 
CALL BT SECUIBDS 

CALL BLOFF(ID2 

XM = XBLK + 0.10 

YM = YBLK + 0.05 

CALL aN 

CALL MESSAG(‘*R(E =)$’,100,XM, YM) 
XR=XM+.5 

CALL REALNO(RE,—1,XR, YM) 


XM = XBLK + @.1 
YM = YBLK + 0.30 

CALL HEIGHT(®. 15) 

CALL MESSAG(’(RF =)$°,100,XM, YM) 
XR = XM + 0.6 

CALL REALNO(RF,3,XR, YM) 


XM = XBLK + 0.1 

YM = YBLK + 0.55 

CALL aN 

CALL MESSAG(’(T =)$’,100,XM, YM) 
XR = XM + 0.6 

CALL REALNO(T1,3, XR, YM) 


XM = XBLK + 0.1 
= YBLK + 0.80 
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oe 








CALL HEIGHT(@. 15) 

CALL MESSAG(’\A) _=$’ , 100, XM, YM) 
XR = XM + 0.6 

CALL REALNO(ALPD,1,XR, YM) 

CALL BLON(ID1) 

CALL BLON(ID2) 


C..PLOT THE CP DISTRISUTION 


CALL WRERY 02) 
CALL MARKER(15) 
CALL CURVE(XU,CPU, IU,@) 
CALL MARKER(18) 


C CALL DASH 
c CALL BE eRe Chen IEXP) 16,8) 
® CALL RESET(’DASH’ ) 
fe CALL CURVE(XL,CPL,IL,1) 
CALL RESET(*THKCRV’ ) 
C CALL GRID(2,1) 
CALL ENDPL(@) 
GO TO 1020 
20@@ ENDIF 


C se#eaeeeee2eE—END OF CP LOOP, BEGIN LOADS LOOPeseceeees 
IF(IOPT.GE.1) THEN 
C..INPUT FROM LOADS 


I=1 
2100 READ(1@,END=22@0) NT1,NTPL,1T1,ALPD,RF,RE,CL1,CD1,CMi,CP 
elt NTs} GO TO 2100 
IF(NT1.GT.NTL) GO TO 2200 
WRITE(6,11) NT1,T1,ALPD,CL1,CD1,CM1 
11 FORMAT(’ °,15,5F8.3) 
ANG( I )=ALPD 
T(1)=T1 
CL(1)=CL1 
CD(1)=CD1 
CM(I )=—CM1 
l=I+1 
GO TO 2100 
2206 CONTINUE 
IPTS=I-1 


C..INNER LOOP — SUS8PLOTS 

IFCIOPT.EQ.1) CALL PHYSOR(3.@,7.75) 
aeer Eo: 22 CALL OREL(@. ,-3.50) 
IF(IOPT.EQ.3) CALL OREL(@. ,-3.5@) 
C..SET UP ALPHA—NUMERICS AND AXES NAMES 

CALL NOBRDR 

CALL COMPLX 

CALL BASALF(’STAND’ ) 

CALL MIXALF(’L/CSTD’ ) 

CALL omescc, 

CALL BLOWUP (SCALE) 
C..LABEL AND DRAW AXES 


IF(I1OPT.EQ.1) THEN 


JE, 


CALL YNAME(’L(IFT) C(OEFFICIENT)$’ , 100) 
Yi=—.5 
N2= 2 
DY=.5 
ENDIF 
IF(IOPT.£Q.2) THEN 
CALL YNAME(°O(RAG) C(OEFFICIENT)$° , 100) 
Yiz=.2 
Y2=1. 
DY=.2 
ENDIF 
IF(IOPT.£Q.3) THEN 
CALL XNAME(’\A)$’ , 100) 
CALL YNAME(’*M(OMENT) C(OEFFICIENT)$’ , 140) 
Yi=—.15 
Y2=.15 
DY=.05 
ENDIF 


CALL AREA20(XLEN, YLEN) 

CALL YAXANG(@.@) 

CALL XTICKS(5) 

CALL YTICKS(5) 

CALL GRAF(X1,0X,X2,Y1,DY, Y2) 


C..PARAMETERS 


IF(IPAR.EQ.1 .AND. IOPT.EQ.1) THEN 
WB = 1.55 
HB = @.55 
XBLK = @.65 
YBLK = YLEN-0.6 
CALL BLTREC(XBLK,YBLK,WB,HB,@.@,4.015) 
CALL BLKEY(1ID3) 
CALL BLOFF(ID3) 
XM = XBLK + @.10 
YM = YBLK + 6.05 
CALL pee eee) 
CALL MESSAG(’R(E =)$’ ,100,XM, YM) 
XR=XM+. 5 
CALL REALNO(RE,—1,XR, YM) 


XM = XBLK + @.1 
YM = YBLK + 0.30 
CALL Beoue (2) 
CALL MESSAG(’(RF =)$°,100, XM, YM) 
XR = XM + @.6 
CALL REALNO(RF,3,XR, YM) 
CALL BLON(ID3) 
ENDIF 


C..PLOT THE LOADS DISTRIBUTION 


QAQQaAQ 


CALL THKCRV(.017) 

IF(IOPT.EQ.1) CALL CURVE(ANG,CL, IPTS,®) 
IF(IOPT.£Q.2) CALL CURVE(ANG,CD, IPTS,@) 
IF(IOPT.£Q.3) CALL CURVE(ANG,CM, IPTS,0) 
CALL DASH 
IF(IOPT.£Q.1) CALL CURVE(ALPHA,CLEXP, 20,2 
IF(IOPT.£Q.2) CALL BU VE MMBRAICS EXOOO eos 
IF(IOPT.£Q.3) CALL CURVE(ALPHA,CMEXP, 2@,0) 
CALL RESET(’*DASH’ ) 

CALL RESET(’ THKCRV’ ) 

CALL THKCRV(.207) 


TS 


C CALL GRID(1,1) 
CALL RESET(' THKCRV’ ) 
CALL THKCRV(.0001) 

€ CALL GRID(2,2) 
CALL RESET(’ THKCRV’ ) 
CALL THKFRM(.@23) 
CALL FRAME 
CALL RESET(’ THKFRM' ) 


C..FINISH SUBPLOTS 


CALL ENDGR(@) 
1OPT=IOPT+1 
IF(IOPT.LE.3) GO TO 2200 
CALL ENDPL(0) 

ENDIF 


C eeseseeenenEND OF LOADS LOOPeeces aces 


oO 


. PLOTTING FINISHED 


3000 CALL DONEPL 
STOP 
END 


Hes 


MO A AANQNANAIANARQAAANANANAAAANAN OO 


PROGRAM CARPET 


PURPOSE: Creates a single plot with multiple Cp vs X/C curves ina 
carpet plot format. The Cp and X/C scales are marked on the 
first curve drawn. The other curves are drawn to the same 
scale, but axes are not drawn although the Cp=0.0 level is 
shown with a dotted line. 


EXTERNAL REFERENCES: 


MODULE DESCRIPTION 

DIP Initializes file im NASA-Ames Device Independent Plot format. 
DISSPLA Graphics software. 

OPENER Prompts for the name of a sequentia! file and opens it. (PROGTOOLS) 
PLFREE Writes "free values” with headings ona plot. fete cis 
PLTYTL Writes a centered character string on the plot. (GRAPHLIB) 


DEVELOPMENT HISTORY: 
DATE INITIALS DESCRIPTION 
Nov. 1986 RCL Original design and implementation. 


AUTHOR: Rosalie C. Lefkowitz Sterling Software, Palo Alto, CA 





IMPLICIT NONE 


¢ Parameter constants: 


INTEGER 

> LUNCRT, LUNDAT, LUNDIP, LUNIN, LUNKBD, MAXI, MAXI2, MAXJ, 

> MXSTR, MAXTIT, NFREE 
PARAMETER 

> ( LUNCRT = 6, 

> LUNDAT = 8, 

> LUNDIP = 10, 

> LUNIN = 12, 

> LUNKBD = 5, 

> MAX I = 40, 

> MAXI2 = 46, IMAXI/2 + 1, 

> MAX J = 20, 

> MXSTR = 81, 

> NFREE = 3 ) 


« Variables: 


REAL 
> ALPHA(MAXJ), CPTH(MAXI,MAXJ), FREVAL(NFREE), XTH(MAXI), 
> x2(3), Y2(3) 


DIMENSION WORK (MAXI) 


REAL 
FRELEN, HIGH, HITFRE, HITITL, WIDE, 
XOFSET, XOR, YOFSET, YOR, YINCH, YFREE, YTITLE, 
XMAX, XMIN, XSTEP, 
CPMAX, CPMIN, CPSTEP 


VVVV 


CHARACTER 
> DATFIL#8@, 
>  FRETXT(3)#(10), 


lhe 


> TITLE*(MXSTR) 


INTEGER 
> I, J, NOGFRE(NFREE) 


DATA 
CPMIN,CPSTEP ,CPMAX/ @. ,-20@. ,-20./, 
DATFIL / ‘CARPET.DAT’ /, 
FRELEN/1@.0/, 

HIGH/3.0/, 
HITFRE/O@.09/, 
HITITL/@. 14/, 
NDGFRE/2, 2, 3/, 
WIDE/2.8/, 

XMIN, XSTEP, XMAX/@.,1.,1./, 
K270..0.,1./, 
XOFSET/.@95/, 
XOR/@.@/, 
YFREE/@.5/, 
YOFSET/. 16/ 
YOR/@.0/, 
YTITLE/6.5/ 


VVVVVVVVVVVV VV VV 


* Free text for plot ltabeling and label for aiphoa vaiue fist: 


WRITE( FRETXT(1), 110@ ) ‘AMPLITUDE’ 
WRITE( FRETXT(2), 110@ ) °ST. MEAN’ 
WRITE( FRETXT(3), 110@ ) ’RED. FREQ.’ 


¢ Initialize plotting device: 


CALL DIP ( LUNDIP, ’CARPET.DIP’, 100 ) 


¢ Prompt for and open data file: 


CALL OPENER ( LUNCRT, 
> ‘Enter data file name (default is CARPET.DAT): °, 
> LUNKBD, DATFIL, LUNDAT, ‘OLD’ ) 


* Reod the data file: 


READ (LUNDAT,1000) TITLE 
WRITE(LUNCRT , 1000) TITLE 
READ Means ( XTH(I),1=1,MAXI ) 
READ (LUNDAT,*) ( FREVAL(1),1]=1,NFREE ) 
DO 10 J=1,MAXJ 

READ Nome ALPHA( J) 

READ (LUNDAT,*) ( CPTH(I,J),1=1,MAXI ) 

1@ CONTINUE 


CALL RESET (’ALL') 

CALL SETDEV ( @, @ ) 

CALL HWROT (’MOVIE’) 

CALL HWSCAL (’NONE’ ) 

CALL NOBRDR 

CALL GRACE ( 0.@ ) 

CALL NOCHEK 

CALL PAGE ( 11.0, 8.5 ) 
CALL BASALF ( ‘STANDARD’ ) 
CALL MIXALF ( *L/CGREEK’ ) 


dips No. 


C « Start with title and free text/free vaiues: 


CALL PHYSOR ( XOR, YOR ) 

CALL AREA2D ( 10.0, 8.@ ) 

CALL PLTYTL ( 0.0, FRELEN, YTITLE, HITITL, TITLE ) 

CALL PLFREE ( NFREE, FRETXT, FREVAL, NOGFRE, HITFRE, YFREE, 
> FRELEN ) 
CALL ENDGR(@) 


C »* Then draw Cp curves: 


XOR = 2.0 — XOFSET 
YOR = 1.0 ~— YOFSET 
DO 50@ J =1,MAXJ 


XOR = XOR + XOFSET 
YOR = YOR + YOFSET 
CALL PHYSOR XOR, YOR ) 
CALL AREA2D (WIDE, HIGH) 
CALL CROSS 
IF ( J.EQ.1 ) THEN 
Cis First piot has full labeled axes: 
CALL XNAME </C 
CALL YNAME "=p", 3) 
GEESE 
GC # Later curves are drawn without labeled axes: 
CALL XNAME : ee J) 
CALL YNAME (’ °, @ ) 
END IF 
CALL GRAF (XMIN, XSTEP, XMAX, CPMIN, CPSTEP, CPMAX) 
CALL CURVE ( XTH(1), CPTH(1,J), MAXI2, @ ) 
IF ( J.GT. 1) THEN 
C Draw dotted pseudo—axes: 
Y2(1) = CPTH(1,J) 
vate = @.@ 
Y2(3) = @.2 
CALL DOT 
CALL CURVE ( X2.,Y2, 3. Om) 
CALL RESET ( ‘DOT’ ) 
END IF 
CALL HEIGHT ( HITFRE ) 
Cs How many inches from the origin is Y= @ ? 
YINCH = YPOSN ( 0.0, 0.2 ) 
C Now print the value of ALPHA: 
CALL REALNO ( ALPHA(J), 2, WIDE + @.2, YINCH ) 
C Identify the alpha value list: 
IF (J.£Q.MAXJ ) 
> CALL MESSAG ('’(a) deg.$’, 100, WIDE + 0.2, YINCH + 0.2 ) 
CALL ENDGR (0) 


500 CONTINUE 
C -» Termination: 
CALL ENDPL ( @ ) 
CALL DONEPL 
STOP ° Normal termination. Plot file is in CARPET.DIP. ’ 
1000 FORMAT ey 
110@ FORMAT (A10) 


END 
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NOTES ON THE EMPLOYMENT OF THE WU CODE AT NASA~AMES 





DISCLAIMER 


The following is a site-specific, step-by-step example intended for the 
novice user to employ the Wu code at the NASA=AMES computational facility. 
It makes no attempt to be all encompassing, but rather tries to provide an 
adequate amount of information to shorten the learning period required to 
gain a minimum working knowledge of the facilities. Much of it can be applied 
to the use of other codes at AMES, but it will have only limited application 
at other sites. 


There are ao number of useful publications which can be helpful in 
utilizing the NASA-AMES computational facilities. Included are: 
1. INTRODUCTION TO VAX/VMS AT AMES 
2. CrayVAX DECnet Interface User’s Guide 
3. A GUIDE TO GENERATING MOVIES USING PLOTSD AND GAS 
4. PLOT3D and PLOT3X Version 3.5 (30 for VAX and 3X for Cray) 
Joan Thompson and Rosealie Lefkowitz usually have these publications. 
Their office is in Bldg 227, Rm 102. Both ladies are very helpful and 
generous with their time. 


You will need access to a VAX account and required password. See an 
NPS representative for this information. ‘Jianps’ is an account on RALph 
currently available for NPS use, which is convenient as DISSPLA and PLOT3D 
are installed on RALph. Individual VAX are also referred to as ‘nodes’. 


The following description assumes the user is at a graphics—capable 
terminal (or ‘green screen’). Use for standard terminals wil! be identical 
with the exception of plotting locations available. Remote use does not 
allow for full screen editing. 

First, logon by typing ‘c *##’ , where ’c’is for connect, and ‘ese’, is for 
the first three letters of the node you are logging onto. You will then be 
prompted for your password (pw). Upper or lower case is allowable for almost 
al! commands. 

All regular keyboard entries are submitted to the VAX via the carriage 
return <cr>, with the exception of command line entries which submitted via 
the ’Enter’ key. 

Now that you are logged into the system, you can check to see what you have 
immediate access to by typing ‘dir’. 

If the file needed is in a subdirectory, type ‘sd [.***]’, where ‘ess’ js 
the subdirectory name. 

Type ‘dir’ (or ‘list’ for added file information) again. 

To edit or to look at the file, type *edt fn’ ,where fn is the filename. 
(i.e. this.file;1 or my.dat) The default version (or number after the 
semicolon) is the highest number. 

You are now in the line editing mode. To go to full screen or keypad 
editing, type’c’ and <cr>. 
The keypad is the at the right end of the keyboard. See the "Introduction to 
VAX/VMS at AMES" for more information. AMES will eventually be converting to 
UNIX (which, like "VAX" is a unique operating system) and have some different 
functions and options than are on “DEC”. DEC and VAX are normally used 


lee 


interchangably, but DEC is properly the name of the manufacturer,while VAX 
is the operating system itself. 
After you are finished changing or looking at the file, enter 
*PF1’ and then keypad '7’. This will put you in the command line mode. 
To save any changes and create a new version of the file, type ‘exit’ enter. 
If no changes have been made or you don’t want to save the changes made, type 
‘quit’ and “Enter”. 
To submit a job to the CRAY X/MP-48, type ‘csub’ and the fn. You will then 
be prompted fer your pw. For additional information, see "Cray VAX DECnet 


Interface User’s Guide”. 


glee No 


The following is an example of a series of Cray and VAX runs to employ the 
Wu code and plotting options. 


NOTES WILL BE INSERTED THROUGHOUT. 


NOTES : FOR CRAY JOB CONTROL LANGUAGE (JCL),AN '#*.' WILL COMMENT OUT A LINE. 
ALL CRAY JCL LINES MUST END WITH A PERIOD. 
THE FIRST LINE OF THE JCL MUST BE THE 'JOB=’ LINE. DON’T HAVE A BLANK 
LINE FOR THE FIRST LINE OR THE FIRST SPACE ON THE FIRST LINE! IT WON'T 
WORK. 
INSERT YOUR OWN JOB NAME (JN), USER ID (US) AND USER PASSWORD (UPW). 
YOUR JN MAY BE REPEATED FOR THE NAME AND VAX ADDRESS (VAX OR NODE 
ADDRESS FOR THIS CASE IS FML. OTHERS INCLUDE RALph,MARs, TOM, etc.) 
*X=’ IS FOR YOUR AMES PHONE EXTENSION. 
AS OF 19APRIL88, THE FIRST TWO TIME LIMIT SIZES ARE T=15@ (CPU SECONDS) 
AND T1800. THE FIRST WILL BE ADEQUATE FOR ALL JOBS UP TO ABOUT 
106 TIME STEPS. THE CPU LENGTH OF EACH TIME STEP IS DIRECTLY PRORTIONAL 
THE SIZE OF THE INVISCID REGION BEING CALCULATED. THEREFORE, PRE=-STALL 
ANGLE OF ATTACK CALCULATIONS WILL BE CLOSER TO .5 CPU SECONDS PER TIME 
STEP AND FULLY STALLED CONDITIONS WILL REQUIRE 1 TO 1.3 CPU SECONDS 
PER TIME STEP. 


THE FOLLOWING IS A POSSIBLE SERIES OF COMPUTER RUNS. 
RUN 1: SAVE DATA FOR STEADY STATE CASE 


JOB , UN=eeeee , T=150. 

ACCOUNT , ACwHenees USzeeees UPWaeeess NAME=eeene X=4269 FML: i seess 

+. 

ACCESS , DN=SENDVAX, 1D=STTROM, OWN=sRFTROM. 

*. 

*. Compile, load and run 

+, 

« .CFT ,ON=CSTAX. (USE THESE OPTIONS FOR ADDITIONAL DEBUGGING INFORMATION) 

CFT ,ON@=A,OFF=CST. (USE THESE OPTIONS FOR ALL REGULAR RUNS. THEY PROVIDE A 

LOR. MINIMUM OF EXTRA OUTPUT INFORMATION. ) 

cre 

NOTE: ENTER YOUR DESIRED DIRECTORY FOR THE OUTPUT DATA TO BE SENT TO. 
SENDING THEM TO SCRATCH WILL HELP AVOID EXCEEDING YOUR ALOTTED 
DISC QUOTA SPACE. YOUR JN CAN BE USED FOR THE ##« BELOW. 


SENDVAX ,ON=FT@1 , VON=’ DUA1: | SCRATCH. ##* | TAPE1.DAT’. (UNCOMMENT THESE 
SENDVAX ,ON=FT@3 , VDN=’ DUA1: | SCRATCH. «** JTAPE3.DAT’. LINES WHEN USING 
SENDVAX , ON=FT 14, VON=’DUA1 : | SCRATCH. e#* | TAPE14.DAT’. THE DISSPLA ROUTINES) 


NOTES: DISSPLA ROUTINES INCLUDE PLOT1,PLOT2 AND PLOTS. 
UNCOMMENT THE FOLLOWING LINE WHEN GENERATING DATA FOR PLOTS3D. 


SENDVAX , ON=FT1@, VON=’ FAR@:[ .NAME]GRID.DAT’. 


*, 
REWIND ,DN=FTQ2. 
SKIPF ,ON=$IN.(BYPASSES ANY SPURIOUS LINES AFTER DATA.) 


NOTE: THE PON IS FOR THE DATA SAVED FOR EITHER THE INITIAL STEADY STATE CASE 
(WHICH IS NEEDED FOR ALL DYNAMIC CASES) OR FOR RESTARTS. 


ACCESS ,ON=FTO7 ,PON=TEST1. 


NOTE: THIS IS THE START OF THE JCL FOR THE SECOND JOB TO BE RUN. SUBMITTING 
MULTIPLE JOBS IS IN THIS MANNER IS CALLED ’JOB CHAINING. ’ 


apo / 


*« .CFT ,ON=CSTAX. 
CFT ,ON=A, OFF=CST. 
LOR. 


*, 

SAVE, DN=FT@8 ,PDN=TEST1.(USE ONLY WHEN THE DATA GENERATED WILL) 
(BE NEEDED FOR A SUBSEQUENT RUN. ) 
SENDVAX , DN=FTO9,, VDN=’DUA1: [SCRATCH. ##* ]TAPEQ.DAT®. 
SENDVAX , DN=FTO4, VDN=’DUA1: [ SCRATCH. ##* ]TAPE4. DAT®. 
SENDVAX , DN=FT11,VDN=’[.***]Q.DAT’. 

/EOF 


@ .15 5.61 (ICST=@ TO COMPUTE THE STEADY STATE CASE.) 
.0001 .0005 .0005 75 1020 
.4 .3 .6 
10 35 5@ 74 45 40 
-1Oe10 25 
20) 2 ece (USE THESE VALUES TO ALLOW A CHECK ON THE OUTPUT SOLUTION. ) 


ICST,RF,ALPS, ICTUR 

WMIN , DFMX , ORMX , KMAX , NCC 
URBI , URBP , URR 
1V1,181,1B2,1V2,NPL,NLB 
DTI ,OTINC ,NTMAX 
NTPL,NTOUT,NTLO 


RUN 2: A SMALL NUMBER OF TIME STEPS TO CONFIRM THE RESULTS ARE AS DESIRED. 
NOTE: ALL THE ABOVE JCL IS THE SAME EXCEPT: 


*. SAVE, DN=FTO8 ,PDN=TEST1.(USE ONLY WHEN THE DATA GENERATED WILL) 
*, (BE NEEDED FOR A SUBSEQUENT RUN. ) 
AND 
* .SENDVAX ,ON=FT10, VDN=’FARO:[ .NAMEJGRID.DAT’. 
4 .15 5.01 (ICST=1,2,3 OR 4 AS DESIRED. ) 
.0001 .0005 .@005 75 100 
.4 .3 .6 
10 35 5@ 74 45 48 


10, 1Ome2o 
29 a 2oeeco (USE THESE VALUES TO ALLOW A CHECK ON THE OUTPUT SOLUTION. ) 


ICST,RF,ALPS, ICTUR 
WMIN, DFMX, DRMX , KMAX , NCC 
URBI , URBP , URR 
IV1,181,1B2,1V2,NPL,NLB 
DTI ,OTINC, NTMAX 
NTPL,NTOUT ,NTLO 


RUN 3: THE FINAL DATA CAN BE SAVED ON TAPE IF A RESTART IS TO BE USED. 


NOTE: ALL THE ABOVE JCL IS THE SAME EXCEPT 


IE eke 


SAVE , DN=FTO8,PDN=TEST1.(USE ONLY WHEN THE DATA GENERATED WILL) 
an (BE NEEDED FOR A SUBSEQUENT RUN. ) 


4 .15 5.01. (ICST=1,2,3 OR 4 AS DESIRED. ) 

.@@@1 .0005 .2005 75 120 

4.3 .6 

10 35 50 74 45 40 

.18 .10@ 58 

5@ 5@ 5@ (SELECT OUTPUT TIME STEP VALUES AS DESIRED. THESE VALUES 
WILL GIVE OUTPUT AND PLOTTING DATA IN THE SAME 4.@ SECOND 
INTERVALS AS THE MANUAL DOES. ALTERNATIVELY, ENTER THE 
SPECIFIC ALLP (OR ALPHA IN DEGREES) VALUES IN ZONST. 
COMMENT OR UNCOMMENT THE LINES ASSOCIATED WITH ALLP 
DEPENDING ON THE MANNER OF CHOOSING THE OUTPUT CONTROL. ) 





ICST,RF,ALPS, ICTUR 

WMIN , OF MX, DRMX , KMAX , NCC 
URBI , URBP , URR 
IVi,1B1,1B2,1V2,NPL,NLB 
OTI,DTINC ,NTMAX 

NTPL ,NTOUT,NTLO 


AFTER EACH RUN, LOOK AT THE OUTPUT FILE TO SEE IF THE JOB HAS RUN TO 
COMPLETION. THIS JS INDICATED IF THE CPU TIME IS REASONABLE FOR THE NUMBER 
OF TIME STEPS RUN, IF THERE ARE NO MESSAGES THAT THE JOB WAS ABORTED, AND 
THAT NORMAL DATA TRANSFER WAS ACCOMPLISHED. THEN, IF IT HAS, PLOTTING CAN 
BE DONE. 


IF THE WU PLOTTING ROUTINES (DISSPLA) ARE TO BE USED ENSURE, THE FOLLOWING 
ARE AVAILABLE IN YOUR DIRECTORY: 


FOR PLOT1: 


1. PLOT1.EXE 
2. PLOT1.DAT 
3. PLOT1.COM 


IF YOU DON'T HAVE A VERSION OF PLOT1.EXE IN YOUR DIRECTORY, COMPLETE 
THE NEXT TWO STEPS. 


COMPILE PLOT1.FOR BY TYPING: 
FOR PLOT? (.FOR AND HIGHEST VERSION 
NUMBER ARE THE DEFAULT.) THIS WILL CREATE PLOT1.OBJ 


LINK PLOT1.OBJ BY TYPING: 
LINK PLOT1,SYS$LIBRARY: INTLIB/LIB,DISSPLA/LIB, INT/LIB 
THIS WILL CREATE PLO1.EXE THIS IS THE FILE NEEDED TO 
ACTUALLY DO THE PLOTTING. 


EDT PLOT1.COM TO ENSURE THAT THE PROPER NAMES ARE INCLUDED IN THE FILE NAMES. 
THESE SHOULD BE THE SAME NAMES AS IN YOUR CRAY JCL. 


YOU SHOULD ALSO HAVE SAVED THE DATA BY UNCOMMENTING THE APPROPRIATE LINES 
IN YOUR JCL. SEE ABOVE. 


PLOT1.COM CAN BE: 

$ GRAPHICS 

$ IF “'*FS$MODE()’" .EQS. “BATCH” THEN SET DEFAULT FARQ:[JIANPS. #=] 
$ DEFINE/USER FOR@Q@1 DUA1:[{SCRATCH. «es ]TAPE1.DAT 
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$ DEFINE/USER FOR@@5 PLOT1.DAT 
$ RUN PLOT1 


FOR PCOT 2: 


1. PLOT2.EXE 
2. PLOTZ.DAT 


3. PL2.COM 


PLOT2.COM CAN BE: 


TO CREATE THE PLOTS AVAILABLE FROM PLOT1, 
OPLOT1 


WHERE **«* IS YOUR JN AND #** IS YOUR DIRECTORY 


TYPE: 


IF YOU DON'T HAVE PLOT2.EXE FOLLOW THE ABOVE DIRECTIONS FOR PLOT?1. 


IF "°° F$MODE()'" .EQS. “BATCH” THEN SET DEFAULT FAR@:[JIANPS. #2] 


ee ¢ 
sen 
aes 
eee 
See 


TAPE . DAT 
TAPES .DAT 
TAPE4 .DAT 
TAPES. DAT 
TAPE14.DAT 


$ GRAPHICS 
$ 
$ DEFINE/USER FOR@@1 DUA1:[SCRATCH. 
$ DEFINE/USER FOR@@3 DUA1:{ SCRATCH. 
$ DEFINE/USER FOR@@4 DUA1:[SCRATCH. 
$ DEFINE/USER FOR@Q@9 DUA1: {| SCRATCH. 
$ DEFINE/USER FOR@14 DUA1:[SCRATCH. 
$ DEFINE/USER FOR@@5 PLOT2.DAT 
$ RUN PLOT2 
$ DIPQMS/DELETE PLOT2.DIP 
TO CREATE THE PLOTS AVAILABLE FROM PLOT2, TYPE: 
OPL2 
FOR PLOTS: 

1. LOADS.DAT 

2. LOADS.EXE 

3. LOADS.COM 

4. PLOT3.EXE 

5. PLOT3.DAT 


LOADS .COM CAN BE: 


PRAHA HH HHH HH 


GRAPHICS 


IF "**F$MODE()’" .EQS. "BATCH" THEN SET DEFAULT FAR@:[JIANPS.e¢. e#] 


DEFINE/USER FOR@@3 DUA1: [SCRATCH. ess 
DEFINE/USER FOR@@4 DUA1:|[SCRATCH. ese 
DEFINE/USER FOR@@5 LOADS.DAT 


RUN LOADS 


| DEFINE/USER FOR@24 NACA.DAT 
DEFINE/USER FOR@1@ FOR@1@.DAT 
DEFINE/USER FOR@14 DUA1: (SCRATCH. «**]TAPE14. DAT 
DEFINE/USER FOR@@5 PLOT3.DAT 


RUN PLOTS 


IF YOU DON’T HAVE PLOTS.EXE OR LOADS.EXE, SEE ABOVE. 


TAPES . DAT 
TAPE4 .DAT 


TO CREATE THE PLOTS AVAILABLE FROM PLOTS, TYPE: 
@LOADS 


200 


FOR USING PLOT3D, THESE ARE A SET OF POSSIBLE DATA ENTRIES: 


GRAPHICS 
PLOT3D 
READ/2D 


THE FOLLOWING ARE RESPONSES TO PROMPTS: 
GRID 


Q 
SUBSETS 
1 88 5 


1 605 


1 


MINMAX =-5 5 =-5 5 
FU 200 

WALLS 

1 8@ 


1 


1 


VECTORS 


P/2D/TEK 


THIS WILL GIVE A REPRESENTATIVE PLOT OF THE VELOCITY VECTORS. BY CHANGING 
THE MINMAX AND SUBSETS VALUES, DIFFERENT ASPECTS OF THE FLOW CAN BE 
HIGHLIGHTED. 


THE LAST COMMAND WILL GENERATE A PLOT ON THE TERMINAL. IF A HARD COPY 
IS DESIRED, AFTER THE PLOT IS VIEWED OR IF YOU ARE SURE YOU KNOW WHAT THE PLOT 
WILL LOOK LIKE, TYPE: P/DIP. DO THIS FOR EACH OF THE PLOTS THAT YOU WISH 
TO SAVE. WHEN FINISHED WITH PLOTSD, EXIT FROM IT AND TYPE: DIPQMS Q. THIS 
WILL QUEVE YOUR PLOTS TO THE LASER PRINTER ON RALPH. USE THE SAME TECHNIQUE 
TO PRODUCE HARD COPIES OF THE DISPLAY PLOTS. THERE IS ALSO A BATCH SUBMITTAL 
FILE FOR PLOT2 AND PLOTS. 
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